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1. Introduction 

The original proposal submitted to NASA under this title entailed a comprehensive 
research program for 36 months duration. In view of the fundamental nature of the proposed 
investigations, NASA elected to support only a revised program for one year and, thus, provide 
the PI with a better starting position for seeking support for the remainder of the program from 
other funding agencies, such as the National Science Foundation. Consequently, the activities 
reported here should be viewed as start-up for ensuing research rather than a rounded program. 


The Revised Statement of Work for the one year effort contained five tasks: 

Task 1. Design and assembly of a state-of-the-art high resolution storage and processing 

system for microscopy images. 

Task 2: Vapor growth kinetics and morphology study with carbon tetrabromide. 

Task 3: Design and assembly of a photothermal deflection setup for vapor concen- 
tration measurements during crystal growth. 

Task 4 : Bridgman growth of iodine single crystals to be used in task 3. 

Task 5: Concentration distribution measurements during iodine crystal growth from the 

vapor. 

2. Work Performed 

In the following we will briefly outline the results obtained under the above tasks. 
Overall, the work performed considerably exceeds the initial projections. In addition, we have 
obtained a breakthrough in the modelling of vapor growth morphologies resulting under 
pronounced anisotropy in attachment kinetics, an activity of great importance for the 
experimental work pursued here, but not planned under the original proposal. 

2.1 Microscopy image storage and processing system 

A system was to be developed for the state-of-the-art storage, digitization, processing, 
display and printing of images from an existing Leitz microscope-interferometer vidicon 
setup.This task has been completed: the system has been designed, assembled and tested. All 
components have at least 1024 x 1024 lines of resolution. In addition to the Dage MTI Hr-2000 
high resolution monitor, that was already available in our Center, the system consists of the 
following components and software packages: 

Compaq Deskpro 286, 640K Ram, 40 Mb Hard Drive, with two 1.2Mb, 1/4" floppy drives 
Compaq VGA high resolution monitor and high resolution Logitech mouse 
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Univision UDC 2600-12M display board 

Tecon DVX 1024 frame grabber 

Dage MTI precision 81 high resolution TV camera 

Pasecon video tube 

Dage 207233-1049/30 sync generator 

Grundig BK224H high resolution video recorder 

Image Pro II Software (IP2UDC) 

Customized software for freeze-frame video to capture high resolution 1024 image. 


The TECON corporation, from which we purchased the majority of the components, has been 
invaluable in the design of this system. They have also assembled and tested the system 
according to our specifications prior to shipment to our Center. 

With this system we are able to video-record (in real time) microscopic images of crystal 
growth morphologies obtained in a growth chamber on the microscope table. These images can 
be recalled in their analog form for the identification of crucial events. Selected images are then 
saved in digital form in the host computer. These data are available for image analysis and 
processing. The overall performance of this system is such that, in combination with the Linnik 
microscope interferometer, a depth resolution of about 50 Angstroms can be obtained. In 
addition, features which are hardly discemable on the original low contrast image, can be 
enhanced for ready recognition. 

Features of the ImagePro software enable us to perform several types of image analysis 
and processing, such as intensity histograms (256 grey scales), line distance and angle 
measurements, spatial frequency filtration as well as contrast changes through stretching and 
sliding of intensity distributions. 

In view of the high resolution of the system and the high cost of high resolution printers, 
we have opted for photographic recording of the (as recorded or processed) monitor images. For 
the cataloging of images we are using HP laserjet prints with 300 dots per inch. 

The total cost of the system (without high resolution monitor) is about $31,000, of which 
$21,000 were provided by this grant. The remainder was provided in about equal parts by NASA 
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grant NAG8-71 1 and by the State of Alabama through the Center for Microgravity and Materials 
Research. 

2.2 Growth kinetics and morphology study with Carbon Tetrabromide 

We have investigated the surface kinetics and morphology of CBr 4 during its growth 
from the vapor by high resolution microscopy and interferometry. Temperatures were chosen to 
include the transition from the low temperature monoclinic to the high temperature cubic phase. 
Rounding of comers and edges as well as the development of new facets has been observed with 
increasing temperature at fixed supersaturation on approaching, but distinctly below the 
transition. The rounding temperatures of crystallographically equivalent comers are equal, and 
different for comers formed by different index faces. As the supersaturation is increased the 
rounded comers and edges resharpen. Such a rounding/sharpening behavior can be alternatively 
induced with successively increasing temperature and supersaturation, respectively. We have 
interpreted this behavior in terms of surface roughening theory and computer simulations (see 
also Monte Carlo modelling below). 

With respect to the polymorphic phase transition of CBr 4 we have found that the high 
temperature phase begins to form at grain boundaries and other macroscopic defects at 
temperatures distinctly below that of the bulk phase transition. The new phase grows partly at 
the expense of the old one, leading to surface depressions around the new growth locations. The 
growth rate of the new phase during the transition increases dramatically with temperature and 
exceeds by far the growth rates observed on the monoclinic and cubic phase alone. These high 
growth rates have been experimentally determined to result from simultaneous deposition from 
the vapor and redistribution of solid material via surface diffusion. This rapid growth is 
associated with morphological instabilities. 

This material will be presented at the 1991 March Meeting of the American Physical 
Society and in a major publication to be submitted to the Journal of Crystal Growth (see 

Sect.2.6). 

2.3 Photothermal deflection vapor growth setup 

The vapor concentration measurement cell developed earlier under NASA Grant NAG1- 
733 was to be improved optically and equiped with a better crystal growth tip. This task has been 
completed. A new design, consisting of an optical quality fluorimeter cell cemented to Pyrex 
glass tubes, was developed. In addition we have improved the flow conditions through the 
upstream iodine reservoir to ensure saturation of the earner gas at low pressure differentials. 
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2.4 Bridgman growth of iodine single crystals 

A melt growth technique was to be developed for single crystals of iodine to be used as 
seeds for the concentration field measurements. This task has been completed. Originally we had 
planned to grow a boule, cleave a seed and mount it to a flat coldfinger. This, however, entails 
considerable handling of the (soft) crystal and thus has the potential for introducing high 
structural defect densities and chemical impurity levels. To circumvent these drawbacks we have 
configured the coldfinger (with a seed selection capillary) as an integral part of the melt growth 
system. After growth of a single crystal from a small volume of melt in and onto the coldfinger, 
it is transferred to the vapor concentration probing setup. 

2.5 Concentration field measurement during growth 

The concentration field about growing iodine crystals was to be determined in three 
dimensions at various vapor temperatures, supersaturations and inert gas flow velocities. 
Concurrent photographic records of the crystal morphology were to be made. This task was 
begun and preliminary, encouraging results were obtained. Under growth conditions, 
concentration measurements can be made with a concentration resolution of 0.2 torr of iodine in 
one atmosphere of nitrogen and a spatial resolution of better than one millimeter in three 
dimensions. At low inert gas flow velocities (twice the diffusion velocity) and high wall 
temperature gradients, we see clear evidence for variations in the concentration distribution due 
to convection rolls. These rolls were also confirmed with flow visualization. At higher flow 
velocities and lower wall temperature gradients the concentration profiles closely resemble 
solutions to the one-dimensional transport equation with the crystal forming an appropriate sink. 

2.6 Monte Carlo Modelling of Growth Morphology with Anisotropic Surface Kinetics 

In addition to the tasks outlined in the Revised Statement of Work (see Sect. 1) we have 
very successfully become involved in the theoretical treatment of the interplay between vapor 
transport (task 5) and anisotropic growth kinetics (task 2) through Monte Carlo modelling. The 
model combines nutrient transport, based on a modified diffusion-limited aggregation process, 
with anisotropic surface kinetics and surface diffusion. Through a systematic variation of the 
simulation parameters (temperature, bond strength and supersaturation), the whole range of 
growth morphologies from fully facetted to side-branched dendritic growth is recovered. 

The diffusion in the bulk nutrient and the anisotropy in the interface kinetics are seen to 
be morphologically destabilizing and stabilizing, respectively. It is found that for a given set of 
simulation parameters and symmetry of the lattice, there is a critical size beyond which a crystal 
cannot retain its stable, macroscopically facetted growth shape. This critical size scales linearly 
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with the mean free path in the vapor. Since both thermal and kinetic roughening reduce the 
kinetic anisotropy, the critical size decreases as either temperature or supersaturation is 
increased. Surface diffusion is seen to stabilize facetted growth on the shorter scale of the mean 
surface diffusion length. 

In simulations with a uniform drift superimposed on the random walk nutrient transport, 
crystal faces oriented towards the drift exhibit enhanced morphological stability in comparison to 
the purely diffusive situation. Rotational drifts with periodic reversal of direction are found to be 
morphologically stabilizing for all faces of the crystal. 

These efforts, which have also been supported by the State of Alabama through the 
Center for Microgravity and Materials Research, have lead to five publications (some invited) 
and several (invited) presentations, which are listed in the following section. 

3. Presentations and Publications of Results obtained under this Grant 

- Rong-Fu Xiao, J. Iwan D. Alexander and Franz Rosenberger, Morphological Evolution of 
Crystals Growing in the Presence of a Uniform Drift: A Monte Carlo Simulation, Phys. Rev. A 
39, 6397-6401 (1989) attached. 

- Rong-Fu Xiao, J. Iwan D. Alexander and Franz Rosenberger, Growth Morhology with 
Anisotropic Surface Kinetics, J. Crystal Growth 100, 313-329 (1990) attached. 

- R.-F. Xiao, J.I.D. Alexander and F. Rosenberger, Simulation of Surface Morphologies in 
Crystal Growth from Vapor, J. Crystal Growth (in print) attached. 

- R.-F. Xiao and J.I.D. Alexander and F. Rosenberger, Growth Morphology of Crystal-Vapor 
Interfaces, Phys. Rev. (in print) attached. 

- F. Rosenberger, R. -F. Xiao and J.I.D. Alexander, Morphological Stability of Interfaces with 
Strong Anisotropy in Growth Kinetics, in Lectures on Crystal Growth, ed. by H. Komatsu, (in 
print) attached. 

- F. Rosenberger, Morphological Stability of Interfaces with Anisotropic Growth Kinetics, 
invited seminar, Forschungszentrum Jiilich, Germany, Division for Solid State Physics, 
September 18, 1990. 
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- F. Rosenberger, Morphological Stability of Interfaces with Anisotropic Growth Kinetics, 
invited seminar. University of Osnabriick, Osnabriick, Germany, Department of Physics, 
September 21, 1990. 

- Rong-fu Xiao, Simulation of 2D Crystal Surfaces Growing from the Vapor, Bull. Am Phys. 
Soc. 35, 308 (1990). 

- Rong-fu Xiao, Simulation of 3D Crystal Surfaces during Growth, Bull. Am Phys. Soc. 35, 
684 (1990). 

- Rong-fu Xiao, Simulation of Surfaces Morphologies in Crystal Growth from the Vapor, 
invited paper at Eighth American Conference on Crystal Growth, Vail, CO, July 1990. 

- Rong-fu Xiao, Simulation of Crystal Growth Morphologies, poster at the Gordon Research 
Conference on Fractals, Plymouth State College, NN, August 1990. 
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Morphological evolution of crystals growing in the presence of a uniform drift: 

A Monte Carlo simulation 

Rong-Fu Xiao, J. Iwan D. Alexander, and Franz 'Rosenberger 
, Center for Microgravity and Materials Research , University of Alabama in Huntsville , 

Huntsville, Alabama 35899 
(Received 20 December 1988) 

A Monte Carlo model is used to simulate the morphological evolution of crystals growing under 
the influence of a uniform drift in the nutrient phase. The model combines nutrient transport (via a 
biased random walk) with anisotropic surface attachment kinetics and surface diffusion. It is found 


that the crystal morphology closely depends upon 
strength), and supersaturation. Facets facing the 
stability as compared with the no-drift situation. 

The growth morphology and morphological stability of 
a crystal is an extremely complex problem, involving, in 
general, an interaction between the effects of nutrient 
transport and surface kinetics. 1-7 Recently, the collective 
effect of nutrient diffusion, based on a modified diffusion- 
limited-aggregation 8,9 (DLA) model, and surface kinetics 
on morphological evolution has been investigated by a 
Monte Carlo (MC) simulation. 10,11 A variety of condi- 
tions, ranging from kinetic- to diffusion-controlled 
growth were examined, and successive transitions from 
compact faceted (surface kinetics limited case) to open 
dendritic morphologies (nutrient diffusion limited case) 
were obtained. 

However, in reality, a condition of pure diffusion in the 
nutrient is rarely found. The trajectories of particles in 
the system are inevitably perturbed in some preferential 
direction in the presence of some external field. Such 
bulk flow may, for instance, be due to buoyancy or sur- 
face tension gradients. There is considerable evi- 
dence 12,13 for the influence of a fluid flow on growth mor- 
phology. In particular, the morphological stability of the 
interface may be altered by the presence of a fluid flow. 

Although progress has been made, the theoretical 
study of the influence of a flow on the morphological sta- 
bility has been limited by the difficulties associated with 
the highly nonlinear governing equations coupled with 
moving boundary conditions. So far, only a few cases of 
crystal growth from the melt with simple boundary con- 
ditions have been solved. 6,14-21 For instance, Delves 6 has 
shown that a forced flow parallel to the interface helps 
stablize the interface against a Mullins-Sekerka instabili- 
ty. 3 In addition, it has been shown that coupling between 
morphological and convective (thermosolutal) modes of 
instability can occur in the presence of nonparallel 
flow. 18-21 

For more complex situations, a microscopic approach, 
e.g., the tracking of individual growth units by MC simu- 
lation 22 appears more promising. In 1983, Meakin, 23 for 
the first time, investigated the drift effect on pattern for- 
mation by applying a MC method. In his approach, a 
drift was superimposed on the random walk of a DLA 
model. A crossover was found from the DLA fractal on 
short length scales to a uniform structure on larger 
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the imposed drift, growth temperature (or bond 
drift direction exhibit enhanced morphological 


length scales as a drift was applied. However, lacking 
realistic elements of attachment kinetics, Meakin’s model 
cannot give insight for a crystal growth system. In addi- 
tion, his computational algorithm cannot readily be 
shown to correspond to a known formulation of macro- 
scopic transport. 

In this communication we use a MC model, based on 
our previous work, 10,11 to study the morphological 
change of a crystal that grows in response to a chemical 
potential gradient. We examine the effect of a nutrient 
drift, and also account for nutrient diffusion and surface 
kinetics. We use a uniform drift velocity superimposed 
on a random walk, in contrast to real systems where the 
bulk velocity depends on the local (viscous) flow 
configuration. Consequently, the result of our formula- 
tion, though providing physical insight, can be con- 
sidered at best qualitative. 

We consider a crystal growth system which consists of 
two basic processes, e.g., nutrient transport and surface 
kinetics. We assume that the motion and aggregation of 
growth units take place on a 2D square grid, and restrict 
our study to a physical system with the following proper- 
ties. (1) The nutrient phase is gaseous. (2) The nutrient 
phase consists of two components; component A, the 
growth species, is highly diluted in an inert gas B , which 
randomizes the motion of A. (3) The system is iso- 
thermal, i.e., heat conduction is rapid compared to mass 
diffusion, thus the latent heat released during crystalliza- 
tion can be ignored. (4) The chemical potential is taken 
to be a linear function of the growth species concentra- 
tion. This results in a Fickian transport equation and 
thus, a random walk may be used to describe diffusion in 
the nutrient phase. 10,11,24 (5) There is a macroscopic as 
well as microscopic drift acting on the individual growth 
unit in the nutrient. 

Linder purely diffusive conditions the nutrient trans- 
port can be described by a simple random walk with an 
isotropic jumping probability. 10,11 However, if there ex- 
ists an externally imposed drift, the individual jumps of 
the random walk will be biased. This is schematically de- 
picted in Fig. 1, where the dashed circle is the source 
from which growth units are released. The vector a is an 
unbiased random walk with equal jumping probability to 

6397 ©1989 The American Physical Society 
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pied nearest neighbors of site i, and f is the interaction 
energy of a molecule with a nearest neighbor. 

If the particle does not stick onto the site it will jump 
to one of its unoccupied neighbor sites. The jump proba- 
bility from site / to a neighboring unoccupied site j (on 
the surface or in the nutrient) is 


r < -i 



( 8 ) 


where c' is the number of unoccupied nearest neighbors 
of site i. Clearly, larger n, results in a higher probability 
that a molecule will jump to site j on the interface. 

The combination of Eqs. (7) and (8) contains the essen- 
tial physics of the surface kinetics for a crystal growing 
from a vapor. It should be realized that the crysta an- 
isotropy has already been taken into account implicitly in 
Eqs. (7) and (8), since n is an anisotropic factor. 

Based on the preceding formulation, we carried out the 
simulation by following the previous procedure including 

the multiple registration technique to reduce the 

noise 10>H In this simulation, a larger size of growth pat- 
tern (10 4 particles) with stable morphology was obtained 
by increasing the mean free path to 10 lattice units. T e 
source boundary was chosen at r s — r max , w ere r ax 
is the radius of gyration of the growing cluster. 1 e 
simulation was implemented on a CRAY X-MP/2 an 
typical cpu times used for each simulation were around 
2000-5000 sec. To show the time sequence of the growth 
patterns, we divided the total growth particles into four 
groups of 2500 each and indicated the morphology asso- 
ciated with each group by a line. 

For comparison, results without a superimposed d 
(“diffusion only”) are presented in Fig. 2. This figure 
shows the effect of bond strength or growth temperature, 
i e <t>/kT on the growth pattern for a constant supersa- 
turation (y = 2.0 or Ajr/kT =0.69). One can see that 
with increasing temperature (or decreasing bond 


strength) the crystal morphology changes. At low tem- 
perature, despite an appreciable supersaturation, the 
crystal is still able to retain a compact faceted form [Fig. 
2(a)l. At higher P’s a branched pattern with four main 
protrusions results [Figs. 2(b) and 2(c)], Finally, at the 
highest value of P, it becomes dendritic. The f ° rmatl ° 
of 8 these different growth morphologies has been ex- 
plained by consideration of the collective effec 
trient diffusion and surface kinetics. 

Figure 3 shows the effect of a drift on growth morpho- 
logies at a fixed drift magnitude (6-50 lattice units), and 
otherwise the same conditions as Fig. 2. One can see t 
as a drift is applied, asymmetric growth patternsappear. 
The sites facing the drift grow much faster than those op- 
posite to the drift. Although the drift affects all cases, 
the response of growth morphology increases with de- 
creasing bond strength or increasing temperature [F 
3(a)— 3(d)]. In addition to modifying the growth rate, the 
drift is also seen to stabilize the faceted growth ,n com- 
parison to the “no-drift” case. Whereas, e.g the no-dnft 
patterns of Figs. 2(b) and 2(c) show deep depress ons in 
the middle of the facets, the corresponding drift cases 
[Figs 3(b) and 3(c)] exhibit much smoother facet surfaces 
on which the flow impinges. The flow causes a stronger 
increase of the mass flux to the center of a facet, as com- 
nared to the edges. This can even lead to the overgrowth 
71 Clonfcee ,h. “hole" in Hg. 3®l ns 

open in the corresponding no-dnft case t^g. ^d)]. 

Figure 4 shows the effect of changes in dnft d recUon 
[Figs 4(a) and 4(b), 5 = 50] and magnitude [Fig. 4(c), 
Fin 4(d) 5=16], One can see that the growth 
pattern depends closely on the characteristics of the e drift. 
As the drift direction changes the crystal changes its 
shape accordingly. The morphologically stabilizing effect 
of the drift on the faceted growth morphology ts^agam 
well illustrated in the sequence a,c,d of Fig. - ’ 

here the supersaturation and the bond strength an /or 
temperature correspond to those of the no-dnft case of 



FIG. 2. No-drift case (K«0). Effect of bond s^ngth and 
temperature on growth pattern with ty/kT-OM. ■ W 

^/kT = 4.60, (b) 0 /kT = 3.9O, (c) ^/kT-2.30, (d) $/kT 
=0.69. 


FIG 3. Effect of drift on growth morphology with 8 50, 

and conditions otherwise the same as Fig. 2. The arrow in i- 
cates the direction of the drift. 
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GROWTH MORPHOLOGY WITH ANISOTROPIC SURFACE KINETICS 

Rong-Fu XIAO, J. Iwan D. ALEXANDER and Franz ROSENBERGER 

Center for Microgravity and Materials Research , University of Alabama in Huntsville , Huntsville, Alabama 35899, USA 
Received 15 August 1989; manuscript received in final form 20 November 1989 


A Monte Carlo model is used to simulate the morphological evolution of crystals growing from an incongruent vapor phase. The 
model combines nutrient transport, based on a modified diffusion-limited aggregation process, with anisotropic surface kinetics and 
surface diffusion. Through a systematic variation of the simulation parameters (temperature, bond strength and supersaturation), the 
whole range of growth morphologies from fully facetted to side-branched dendritic growth is recovered. The diffusion in the bulk 
nutrient and the anisotropy in the interface kinetics are seen to be morphologically destabilizing and stabilizing, respectively. It is 
found that for a given set of simulation parameters and symmetry of the lattice, there is a critical size beyond which a crystal cannot 
retain its stable, macroscopically facetted growth shape. This critical size scales linearly with the mean free path in the vapor. Since 
both thermal and kinetic roughening reduce the kinetic anisotropy, the critical size decreases as either temperature of supersaturation 
is increased. Surface diffusion is seen to stabilize facetted growth on the shorter scale of the mean surface diffusion length. In 
simulations with a uniform drift superimposed on the random walk nutrient transport, crystal faces oriented towards the drift exhibit 
enhanced morphological stability in comparison to the purely diffusive situation. Rotational drifts with periodic reversal of direction 
are found to be morphologically stabilizing for all faces of the crystal. 


1. Introduction 

The characterization of the conditions under 
which a growing crystal is capable of preserving 
its shape, i.e. is morphologically stable, is both 
scientifically challenging and technologically im- 
portant. Morphological stability is necessary for 
the growth of homogeneous single crystals that are 
needed for numerous device applications. To ade- 
quately describe the growth morphology and mor- 
phological stability of a crystal, however, is an 
extremely complex problem, involving, in general, 
an intraction between nutrient transport and in- 
terface kinetics [1-3]. 

Two basic approaches to the morphological 
description of crystal growth have been adopted. 
The first is the macroscopic approach which in- 
volves the solution of a continuum transport equa- 
tion coupled with moving boundary conditions [4]. 
The second is the simulation of microscopic 
processes by tracking the individual growth units 

[5]. 

The continuum approach has lead to significant 
insight into the morphological evolution of essen- 

0022-0248/90/503.50 © Elsevier Science Publishers B.V. 
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tially isotropically responding, i.e. atomically 
rough and, hence, macroscopically non- facetted 
[3] interfaces, such as prevail in many melt growth 
systems. Morphological stability conditions are 
well established for a variety of simple geometries 
[4,6-9], including systems with weakly nonlinear 
[9-12] and highly nonlinear isotropic response 
[13-15]. In addition, the coupling between mor- 
phological and hydrodynamic instabilities has been 
examined [16-22]. Continuum models have also 
been developed for the morphological stability of 
non-facetted interfaces with anisotropies in surface 
tension and growth kinetics (e.g., refs. [23-26]). 
One of the most important results of these treat- 
ments is the insight that anisotropies can stabilize 
otherwise unstable closed growth forms up to a 
certain critical size. For atomically smooth inter- 
faces, however, that prevail in most vapor and 
solution growth systems, and on which strongly 
anisotropic (facetted) growth occurs via atomic 
layer spreading, the necessary conditions for mor- 
phological stability are only partly understood. 
Several workers have proposed an isotropic con- 
tinuum formulation for morphological stability in 
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interface via diffusion and convection in the nutri- 
ent. The actual transport kinetics will, in general, 
be determined by the interaction of growth units 
(or their precursors) and other species that form 
the nutrient. Second, when a growth unit reaches 
the interface, typically it does not become im- 
mediately incorporated into the growing crystal. It 
will adsorb and diffuse on the interface in an 
attempt to find an energetically favorable “final” 
attachment site, or it will even return to the nutri- 
ent before it finds such a site. The latter happens, 
for instance, when the growth unit is misoriented 
and cannot form stable bonds with the crystal. 
Alternatively, the growth unit may impinge onto a 
site with too few neighbors to prevent it from 
being dislodged by thermal vibrations before it 
becomes adsorbed. Or, a growth unit may even 
become dislodged after having arrived at an en- 
ergetically favorable (“final”) site. The probability 
of all these steps is determined by the local config- 
uration (number of bond-forming neighbors) of 
the interface sites that the unit happens to visit 
during its (short) residence time. At molecular 
length scales, the crystal symmetry will inevitably 
be reflected in the anisotropy of the attachment 
kinetics. This underlying symmetry will also be 
manifested at macroscopic length scales unless the 
atomic roughness of the crystal interface is high 
enough to allow diffusion in the nutrient phase to 
control the evolution of the growing shape. Though 
a complete model of such complex scenarios is not 
practical at this point, we will formulate a Monte 
Carlo model which retains the essential physics of 
both nutrient transport and interface kinetics, in- 
cluding surface diffusion. 

2.1. Mass transport in the nutrient 

We assume that the gaseous nutrient phase 
consists of two components: a growth species A, 
highly diluted in an inert gas B, such that the B 
concentration is essentially uniform and A-A in- 
teractions can be ignored. Component B rando- 
mizes the motion of A, though both components 
may be subjected to a uniform drift. In addition 
we assume that heat dissipation is rapid compared 
to mass diffusion; thus the latent heat released 
during crystallization can be ignored and isother- 



Fig. 1. Schematics of 2D random walk at a large mean free 
path. Dashed circle: source. Dotted central region: growing 
crystal. Concentric dotted circles: contain possible landing 
points (with automatic rounding, i.e. |a|±0.56). The two 
small collision zones have been enlarged to illustrate the differ- 
ent collision behavior in the bulk and the interface. 

mal conditions prevail. The chemical potential is 
taken to be a linear function of the growth species 
concentration. This results in a Fickian transport 
equation and, thus, a random walk may be used to 
describe diffusion in the nutrient phase [47-60]. 

These assumptions are consistent with the fol- 
lowing Monte Carlo model. The motion of a given 
growth unit in the nutrient is described by a 
simple random walk on a discrete lattice of spac- 
ing b, with equal jump length, \a\, in all direc- 
tions, unless a (uniform) drift is superimposed on 
the random walk. This is schematically depicted in 
fig. 1. The large dashed circle of radius r s is the 
source from which walkers are released. After 
release, a random walker jumps with equal prob- 
ability to all grid points within the annular region 
given by | a | ±0.5Z> (with automatic rounding to a 
nearest grid point if the chosen jumping site is not 
exactly on a grid). The jump length | a | represents 
a physical distance on the order of the mean free 
path of molecules. The grid spacing b y on the 
other hand, represents a lattice unit of the crystal 
which is denoted by the dotted central area. In a 
real crystal-vapor system the mean free path of a 
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a random walk process, the chemical potential of 
the source does not appear explicitly in our model. 
The Ap in eq. (4), therefore, is governed by the 
difference between the bulk transport-dependent 
vapor concentration at the interface and the equi- 
librium concentration at the same temperature. 
For a more detailed discussion, see ref. [58]. 

After impingement on the surface, growth units 
can either remain at the original site, return to the 
vapor (evaporate) or continue to wander to an 
unoccupied neighboring site (surface diffusion). 
The rates of both the evaporation and surface 
diffusion processes are sensitive to the local con- 
figuration of the site from which a unit is to be 
dislodged. Hence, following Gilmer and Bennema 
[46], we cast the evaporation rate into the site-de- 
pendent form 

K~ = v exp( - E,/kT ) , (5) 

where v is a lattice vibration factor and E t is the 
product of the pair interaction (bond) energy <j> of 
a unit with a nearest neighbor and the number of 
occupied neighbor sites of site i, n r 

The probability that a growth unit sticks onto a 
specific site is then P, = K' /(K' + K, ). In order 
to relate and v we assume local equilibrium, 
and equate impingement and evaporation rates at 
equilibrium and obtain [58] 

P. = + yP n «-”‘), ( 6 ) 

where w 0 , the number of nearest neighbors in a 
kink site, is 2 for a square lattice and 3 for a 
triangular lattice, and y = exp( Ap/kT) and /? = 
exp( — <f >/kT). 

If the growth unit does not stick onto the site 1 , 
it will jump to one of its unoccupied neighbor 
sites. The jump probability from site i to a 
neighboring unoccupied site j (on the surface or 
in the nutrient) is assumed to be [58-60] 



(7) 


where c f is the number of unoccupied nearest 
neighbor sites of site /. Clearly, a larger rij results 
in a higher probability that a molecule will jump 
to site j on the interface. 

The combination of eqs. (6) and (7) contains 


the essential physics of the surface kinetics for a 
crystal growing from a vapor. It should be realized 
that the crystal anisotropy has already been taken 
into account implicitly in eqs. (6) and (7), since 
both n i and rij are anisotropic. 

The strength of surface diffusion (described by 
eq. (7)) is characterized by the diffusion time or 
diffusion length on the surface. An exact calcula- 
tion of surface diffusion length is difficult. In our 
model we estimate the average surface diffusion 
length X s according to 

where n t is the average (residence) “time” (actu- 
ally the number of MC steps) for a random walker 
to diffuse on the interface (see section 2.3). 

2.3. Simulation procedure 

A summary of the steps involved in the MC 
simulation is presented in fig. 2. Initially, a small 



Fig. 2. Flowchart of simulation steps. 
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Fig. 3. Effect of bond strength and temperature on growth 
patterns of a crystal (triangular lattice) at a fixed (normalized) 
supersaturation (Ap/kT= 0.69): (a) <p/kT= 3.91; (b) $/kT 
= 2.30; (c) 4>/kT= 0.69; (d) <}>/kT= 0.36. 





lines) around growing crystals, (a) circular crystal; (b), (c), (d) 
growth conditions corresponding to the terminal sizes of figs. 
4a, 4c and 4d, respectively. 


under these or even lower <f>/kT conditions, leads 
to the evolution of “macroscopic” depressions 
(figs. 3b and 4b). On further increase of the tem- 
perature, this transition occurs at even smaller 



Fig. 4. Effect of bond strength and temperature on growth 
patterns of a crystal (square lattice) at a fixed supersaturation 
{Ap/kT= 0.69): (a) <f>/kT = 4.60; (b) 4>/kT= 3.91; (c) 
*/kT= 2.30; (d) 4/kT = 0.69. 


sizes and dendritic growth occurs (figs. 3c and 4d) 
which subsequently exhibits extensive side-branch- 
ing (fig. 3d). 

The concentration distributions about a cir- 
cular growth pattern and the terminal crystal 
shapes of figs. 4a, 4c and 4d are presented in figs. 
5 and 6. The iso-concentration lines in fig. 5 are 
spaced by constant concentration increments. 
Hence, more closely spaced lines indicate steeper 
concentration gradients. This is further illustrated 
in fig. 6 by the corresponding concentration pro- 
files along the radius vector through a comer and 
center of a facet, respectively. As expected, for the 
round pattern (test case for our algorithm) the 
concentration distribution is circular, with the 
concentration gradient decreasing with increasing 
distance from the growing interfaces. For the 
“crystals”, however, one sees that both concentra- 
tion and concentration gradient are significantly 
higher in the comer regions than at the center of 
the facet (or what used to be a facet at smaller 
crystal size), in agreement with numerous experi- 
mental findings [62-64]. These differences in- 
crease as the growing crystal takes on a more 
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shape preservation for facetted crystals [36] and its 
adaption by Kuroda et al. [37]; see also ref. [42]. 
The stabilizing effect of the increase in step den- 
sity as the face center depression increases was 
expressed in terms of kinetic coefficients in 

ref. [36] and B( p) in ref [42]), where p is the local 
slope of the vicinal surface segment (local tangent 
to the depression). The growth rate normal to the 
macroscopic face (fictitious singular surface), V n , 
is then related to the nonuniform, normalized 
supersaturation, a, in the form 

V n = /S(p)o. ( 9 ) 

It is interesting to re-examine the microscopic 
morphologies assumed to be relevant in these 
works, which postulate equivalence of the kinetic 
effect of microscopic depressions (fig. 7a) and 
protrusions or hillocks (fig. 7b). Depressions are 
thought to prevail at higher supersaturations, 
where 2D nucleation in the comer regions is likely 
to form the dominant source for growth steps. 
Protrusions, on the other hand, are thought to 
prevail at lower supersaturation, where defects, 
such as screw dislocations (that supposedly form 
more readily in the center region) supply the 
growth steps. These features have been well char- 
acterized for stable growth from a more uniform 
nutrient field. It is not clear, however, how an 
initially shallow hillock, in order to become kineti- 
cally more active, can increase in central height 
when, as growth proceeds, the nutrient supply 
progressively dwindles in the center region. Our 
simulations, as well as numerous growth-mor- 
phology investigations of actual crystals [67], do 
not lead to protrusions at face centers. 

The above simulations also show, for the first 
time, the destabilizing effect of an increase in the 
growth temperature. This is revealed explicitly in 
figs. 3 and 4 which present a transition from 
facetted to non-facetted shapes as the temperature 
is increased. At low temperatures the morphology 
is controlled by surface kinetics while at higher 
temperatures volume diffusion is seen to be 
shape-determining. Fig. 8, which has been com- 
puted from (6), shows that at a fixed supersatura- 
tion an increase in temperature results in reduc- 
tion in the anisotropy of the sticking probability. 



NORMALIZED BOND STRENGTH, */ kT 

Fig. 8. Dependence of sticking probability P t (eq. (6)) on 
normalized bond strength/temperature and number of oc- 
cupied neighbor sites at fixed Ap/kT — 0.69 for (a) square 
lattice (n 0 - 2) and (b) triangular lattice (n 0 = 3). 

The temperature (or bond strength) independence 
of P 2 in fig. 8a (P 3 in fig. 8b) results from the 
assumption of local equilibrium underlying (6), 
where the vapor-crystal equilibrium is with re- 
spect to the kink site. This reduction in the dif- 
ferences between the P/s causes the growth pat- 
tern to become more isotropic, more dependent on 
the concentration distribution, and thus more 
volume diffusion controlled. Note that this transi- 
tion is, loosely speaking, analogous to thermal 
roughening [68,69], The latter, however, is associ- 
ated with equilibrium conditions and is formally 
based on minimization of the crystal’s total surface 
free energy under the constraint that the volume 
be constant. This constraint is not applicable to 
the dynamic growth conditions considered here. 

Furthermore, a comparison of figs. 3 and 4 
shows that a decrease in symmetry of the lattice 
decreases the ability of the crystal to retain a 
facetted growth form. For instance, facetting pre- 
vails at <t>/kT= 3.91 on the triangular lattice (fig. 
3a). Yet, the perimeter of the square crystal shows 
macroscopic depressions at the same normalized 
temperature and supersaturation (fig. 4b); a fully 
facetted form only appears at significantly higher 
values of <f>/kT (fig. 4a). This difference in stabil- 
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0 12 3 4 


NORMALIZED BOND STRENGTH , q>,/kT 

Fig. 10. Crystal-growth morphology as a collective effect of 
supersaturation and bond strength (or temperature) for trian- 
gular lattice, Af = 3xl0 3 and mean free path of one lattice 
unit. In region I crystals acquire compact- facetted forms; in 
region II compact branched with six-fold symmetry; in region 
III dendritic with multiple sidebranches. Symbols: Monte Carlo 
results judged as “boundary cases”. Solid line: Temkin’s 
boundary between facetted and non-facetted growth for a 
simple cubic lattice [78]. For more details, see ref. [58]. 


loss of morphological stability, i.e. under bulk 
diffusion-controlled conditions, any macroscopi- 
cally facetted growth form will become unstable 
beyond a certain M . 

3,2. Changes in growth morphology with increase of 
mean free path 

The above results are all based on the assump- 
tion of a mean free path equal to one lattice unit, 
i.e. | a | = h. In this section we increase the mean 
free path and simulate growth patterns containing 
10 4 particles. To show the time sequence of growth 
shapes, the total particles are divided into four 
groups of 2500 each. 

Fig. 11 shows the effect on morphology of 
bond strength or temperature at a constant nor- 
malized supersaturation, a nd a mean free path of 
11 lattice units with otherwise the same growth 
conditions (i.e., and y) as fig. 4. It can be seen 
that the critical size of a facetted crystal has 
increased dramatically. In fig. 4b the crystal has 
lost its macroscopically facetted form at a size less 
than 3000 particles, but in fig. lib the facetted 
crystal remains up to around 6000 particles. In fig. 
4d the crystal has become a dendrite with many 


son for this is that nutrient diffusion was not 
taken into account by Temkin. As pointed out 
above, nutrient diffusion is always a destabilizing 
factor and reduces the range over which facets can 
exist. 

It should be emphasized, that the specific loca- 
tions of the boundaries between the different mor- 
phology regions presented in fig. 10 are only rep- 
resentative for growth patterns containing M = 3 
X 10 3 particles. For smaller sizes the transitions 
will occur at higher supersaturations, and/or lower 
bond strength (higher temperature). For larger 
sizes, on the other hand, loss of morphological 
stability would set in at combinations of lower 
supersaturations and higher bond strengths than 
in fig. 10. Hence, there is some scaling between 
the overall size of a growth pattern and its critical 
size. Thus, even under most favorable growth con- 
ditions, any system will eventually reach a size 
where the destabilizing bulk diffusion will cause 



(a) 





Fig. 11. Effect of bond strength and temperature on growth 
patterns at a large mean free path ( | a | =10 lattice units). The 
growth conditions (y and /?) are the same as in figs. 4a-4d. 
Contours correspond to addition of 2500 particles each. 
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Fig. 14. Dependence of critical size on mean free path at 
8p/kT~ 0.69: (a) <t>/kT= 4.60; (b) $/kT= 3.91; (c) <f>/kT 
= 2.30; (d) 4>/kT= 0.69. 


simulation of critical crystal sizes determined ex- 
perimentally [31-35]. However, fig. 14 encourages 
us to scale linearly to these sizes. For instance, 
Nanev and Iwanov [33] have found a critical size 
of 120 fx m for a Zn crystal (6 — 1.38 A) growing 
from the vapor phase consisting of a Zn partial 
pressure of approximately 0.1 Torr in 150 Torr of 
hydrogen. The mean free path of H 2 at this pres- 
sure and temperature is about 1.2 jx m. In order to 
obtain a critical size for a relatively large <\>/kT 
system, such as the above Zn experiment, we have 
run fig. 11a beyond the terminal size shown. 
“Macroscopic” loss of face stability was observed 
at a size of 1306, i.e. with b = 1.38 A at 1.8 X 10 -2 
jam. Remember, that fig. 11a is based on \a \ = 
106, i.e. a mean free path of 1.4 x 10" 3 /im was 
assumed. Hence, it is not surprising that the criti- 
cal size of a face obtained from our MC simula- 
tion is three to four orders of magnitude smaller 
than the experimental finding. 

3.3. Surface diffusion and facet stability 

The effect of surface diffusion on growth mor- 
phology has been studied both theoretically 
[1,46,58,81] and experimentally [81]. As might be 
expected, it was found that surface diffusion can 
smooth the growth morphology on the scale of the 
average surface diffusion length A s . This can be 


explained as follows: Wherever the surface has a 
high curvature the nutrient concentration and its 
gradient are also high. This is shown in figs. 5 and 
6. When surface diffusion is permitted, the con- 
centration gradient along the interface is the driv- 
ing force for surface flow from a protuberance to 
a depression. This makes the protuberance less 
stable, while stabilizing the original facet. 

The contribution of surface diffusion to face 
and facet stability is a function of the surface 
diffusion length A s . For the stability of the entire 
face, surface diffusion is important only when A s 
is comparable to the face length. Hence, the effect 
of surface diffusion is particularly evident in our 
simulations when | a ] =6. However, as a crystal 
grows larger, the size of faces eventually exceeds 
A s . Surface diffusion will then continue to play a 
stabilizing role only for terraces or facets with 
sizes less than or on the order of the average 
surface diffusion length. We have found that even 
for the dendritic growth, surface diffusion still has 
a noticeable effect on primary branch thickening 
[58]. 

The above trends can clearly be deduced from 
figs. 3, 4, 11 and 12. For a quantitative evaluation, 
however, one must take into account that the 
surface diffusion length decreases with increasing 
temperature (decreasing bond strength) and super- 
saturation, while being insensitive to changes in 
the mean free path. Consequently, as can be seen 
from figs. 11 and 12, the terrace size decreases in 
response to an increase in temperature, but re- 
mains approximately constant as the mean free 
path is increased. The connection between the 
temperature dependence of A s and face stability is 
further illustrated in fig. 15, which shows a 
dramatic decrease in the average surface diffusion 
length on loss of the planar growth morphology. 
Therefore, particularly at low temperature (high 
bond strength) and low supersaturation, surface 
diffusion must be taken into account in stability 
considerations. 

3.4. Changes in crystal morphology 4n the presence 
of a drift 

With an externally imposed drift, i.e., d ¥= 0 in 
(1), the individual jumps of a random walker will 
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(a) 




(b) 




Fig. 18. Effect of periodically reversed rotational drift (\d\ = 
2b) on growth morphology: (a), (b) Ap/kT— 0.69, 4>/kT = 
3.91; (c), (d) Ap/kT = 0.69, <}>/kT=230. The rotational di- 
rection is reversed after every 2500 particles in 9a) and (c), and 
500 particles in (b) and (d). 


deeply depressed morphologies result. This is in 
agreement with the common observation of the 
benefits of frequent rotation reversal in solution 
growth [83,85,86]. 
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Note added in proof 


ing morphological changes under an influence of a 
uniform nutrient rotation of | d | = 2b at | a | = 
106; the growth conditions are otherwise the same 
as in figs. 11 and 16. The center of rotation 
coincides with that of the crystal. One an see that 
the rotational drift causes the depressions to be 
shifted away from the center of the facet against 
the direction of the drift. This is, for instance, 
observed in solution growth where uni-directional 
rotation causes depressions and “ veiling” (occlu- 
sion of mother liquor) behind the leading edges of 
the crystal; see, e.g., fig. 39 in ref. [83], and ref. 
[84]. 

Fig. 18 shows the morphologies obtained after 
the rotational drift has been periodically reversed. 
The growth conditions (y and ft) of figs. 18a and 
18b correspond to fig. lib and figs. 18c and 18d 
to fig. 11c. The rotation direction is changed after 
the addition of each set of 2500 growth units in 
figs. 18a and 18c, while in figs. 18b and 18d the 
direction is reversed after every 500 particles. One 
can see immediately that crystals with stable mor- 
phology result when to rotation is reversed more 
frequently (figs. 18b and 18d), whereas in the 
corresponding no-drift cases (figs, lib and 11c) 


Most recently, Saito and Ueta have obtained 
similar growth shape evolutions from a Monte 
Carlo model with anisotropic attachment and 
evaporation kinetics features, using a DLA tech- 
nique which launches many walkers at once [87]. 
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ABSTRACT 

The morphological evolution of crystals that grow from an incongruent vapor by the 
surface nucleation and screw dislocation mechanisms is simulated with a Monte Carlo model. 
The model combines volume transport, based on a modified diffusion-limited aggregation 
process, with anisotropic surface kinetics without a SOS restriction. It is found that for a given 
set of simulation parameters and symmetry of the lattice, there is a critical size beyond which a 
crystal cannot retain its stable, macroscopically faceted shape. This critical size scales linearly 
with the mean free path in the vapor. While surface diffusion is seen to stabilize the growth 
morphology on the scale of the surface diffusion length, volume diffusion is always 
destabilizing. Surface roughness increases with increase in growth temperature and 
supersaturation, which reduce the anisotropy in kinetics through thermal and kinetic 
roughening, respectively. For the screw dislocation mechanism, that can dominate at low 
temperature and supersaturation, we find that the combined effect of bulk and surface diffusion 
reduces the terrace width of a growth spiral in its center region. At elevated temperature and 
supersaturation normal growth can dominate in comer and edge regions of a crystal, while the 
spiral growth mode prevails in the center of a facet. 
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1. Introduction 

Monte Carlo simulation has become a popular method in the study of a crystal growth 
process since Chernov and Lewis [1]. Typical studies include equilibrium and growth 
morphologies of crystals [2,3], surface roughening transitions [4,5] and growth rate 
dependence on supersaturation and temperature [6,7]. For reviews, see [8,9]. These studies 
were primarily focused on interfacial kinetics and have occasionally included surface diffusion 
[7]. A solid-on-solid (SOS) restriction [10] (no overhangs) was often assumed. The influence 
of volume diffusion, i.e., the transport of growth units to the interface, has received less 
attention until recently [11-15]. In the event that surface diffusion and interfacial kinetics 
govern the growth morphology, these simplifications are not severe limitation. However, in 
reality, volume transport often plays a decisive role in limiting morphological stability [16-18]. 

We have developed a MC model to study the morphological evolution of growing 
crystals by considering both volume transport and surface kinetics [11-13]. For the volume 
transport, we modified the diffusion-limited aggregation (DLA) model of Witten and Sander 
[19]. Our results showed that the diffusion in the volume and the anisotropy in the interface 
kinetics are morphologically destabilizing and stabilizing, respectively. Through a systematic 
variation of the simulation parameters (temperature, bond strength and supersaturation), the 
whole range of growth morphologies from fully faceted to side-branched dendritic growth was 
recovered. For a given set of simulation parameters and symmetry of the lattice, there is a 
critical size beyond which a crystal cannot retain its stable, macroscopiclly faceted shape. In 
this short paper we will pay special attention to the scaling of the critical size as a function of 
the mean free path in the nutrient. We will investigate the dependence of surface roughness on 
temperature or bond strength and supersaturation. For crystals with a screw dislocation, we 
will demonstrate how surface and volume diffusion influence the terrace width of the spiral. In 
addition, we will study the competition between the two-dimensional nucleation (2DN) and 
screw dislocation growth mechanisms as a function of temperature and supersaturation. 
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2. Model 

We assume the vapor to be isothermal and composed of an inert gas B in which a 
growth species A is highly diluted. In addition, A-A interactions, convection and latent heat 
released during crystallization are ignored. A random walk is used to describe the diffusive 
motion of growth units from the source [11]. The complex process following the arrival of 
growth units at the interface is described by impingement, evaporation and surface diffusion 
rate equations [11-13]. The impingement rate K+ is assumed to be only dependent on 
supersaturation (chemical potential difference A|i) between the nearby nutrient and the crystal 
surface, i.e. K + = K eq exp(A|i / kT) with Keq being the temperature dependent equilibrium 
value of K + . The evaporation rate, on the other hand, is assumed to be sensitive to the local 
configuration of the site from which a unit is to be dislodged. Hence, following Gilmer and 
Bennema [7], we cast the evaporation rate K- into the site-dependent form Kj = v exp(- 
Ej/kT), where v is a lattice vibration factor and, in a nearest neighbor approximation, Ei is 
simply the product of the pair interaction (bond) energy <(> of a unit with a nearest neighbor and 
ni the number of solid neighbors of site i. For the surface diffusion, we assume that the 

diffusion rate depends on the occupation condition of both the site i that the particle occupies 
and the potential jump site j. Hence, we express the jump rate as K^j = v s exp(-AEij /kT). 

Here Vs is a surface vibration factor, and the activation energy AEy = <j>(nj -nj) +8j, where is 
the activation energy when nj ^ nj. With the above rate equations, the corresponding 
probabilities for attachment, evaporation and surface diffusion can be easily obtained [11,13]. 
The details of the simulation procedures have been reported elsewhere [11-13]. 

3. Results 

The results in fig.l were simulated on a two dimensional (2D) square lattice. Initially, a 
nucleus with a certain size was put at the center of a circular periphery from which random 
walkers (growth units A) were released. The radius of this source periphery was chosen such 
that the growth shape was not biased [11]. In this figure, both the temperature and 
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supersaturation were kept constant (<J>/kT = 0.69, A|i/kT = 0.69), and only the mean free path 
|a| was changed respectively from lb, 5b, 10b, 25b, 50b to 100b, with b being the lattice 
constant. Also, a multiple registration scheme [11] was employed to reduce microscopic noise. 
As can be seen from this figure, a crystal can have very different morphologies when the mean 
free path is changed. At |a| = b (fig. la) the crystal shows dendritic patterns with many side- 
branches from the earlier stages of growth. As |a| increases less side-branches are observed 
(fig.lb,c), and when the mean free path increases to 25b (fig. Id) the side branches disappear, 
leaving only depressions at the face center. At a mean free path of 100b, the crystal becomes a 
square without any depression. Similar changes of the crystal morphology in three dimension 
(3D) are shown in fig.2. In this case, instead of a circular source, a spherical source was used 
[13]. Due to computational time limitations, only two different mean free paths were 
considered, i.e., |a| = lb (fig.2a,b) and|a| = 5b (fig.2c,d). According to gas kinetic theory, the 
mean free path is inversely proportional to the total vapor pressure. An increase in the mean 
free path implies a decrease in the vapor pressure or diffusion bamer. Figs.l and 2 
demonstrate the importance of the total vapor pressure for the morphological stability in vapor 
crystal growth. 

Our results also emphasize a fundamental difference between crystal growth and DLA 
[19]. As can be seen from the evolution of the growth boundaries in fig. 1, the scale invariance 
exhibited by DLA does not generally hold in crystal growth. The instantaneous fractal 
dimension decreases as the size of the growing crystal increases. For a given set of growth 
conditions, a small crystal may be stable, but as it grows bigger and exceeds a certain critical 
size, the original crystal shape can lose its stability, and develop a different morphology. This 
observation is compatible with the earlier theories [9,18] of faceted growth shape stability. 

To quantitatively characterize the critical size of a growing crystal at the onset of 
instability, the change of total surface area has been used as a criterion [13]. For a stable 
morphology the total surface area (or total surface length in 2D) should linearly scale as a 2/3 
power (or 1/2 in 2D) with the number of growth particles. When a crystal loses its stability , the 
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surface area increase significantly [13]. Fig.3 shows the dependence of the critical size on the 
mean free path at fixed Ap/kT for both two and three dimensions. As can be seen, the critical 

size increases with decreasing temperature. At fixed temperature and supersaturation, the 
critical size scales linearly with the mean free path in the parameter range considered. The 
corresponding 3D simulations yield larger slopes, which is due to the fact that, on the average, 
there are more solid neighbors associated with interfacial particles in three dimensions and, 
thus, the (stabilizing) anisotropy in interface kinetics is more pronounced. Although 
computational time limitations do not allow for the direct simulation of the critical size at the 
mean free path used in morphological stability experiments [20], fig.3 encourages us to scale 
linearly to these conditions. Order-of-magnitude agreement is obtained between experimental 
and modelling results [12] for the critical size. 

Fig.4 shows the effect of surface diffusion on the morphology of a growing crystal. 
This figure was simulated with a special geometry consisting of two parallel infinite planes, of 
which one is the growing crystal and the other is the source. Periodic boundary conditions 
were used laterally in both x and y directions. In this case, instead of using the multiple 
registration scheme to reduce the noise [11], already attached particles were randomly allowed 
to evaporate at a rate controlled by a Boltzmann factor [13]. Volume diffusion is not considered 
in this figure, i.e. we have taken the ballistic impingement approach used, e.g., by [7]. For 
comparison, a result without surface diffusion is plotted in fig.4a, while fig.4b includes 
surface diffusion under otherwise the same growth conditions (<>/kT=1.6, A[i/kT=0.69). It can 
be seen that with surface diffusion the crystal surface is much smoother (fig.4b). This is 
because surface diffusion can provide an additional way for interfacial particles to relax to 
some energetically more favorable (low energy) sites and, hence, smooth the surface 
morphology. As we have discussed before [12], such a smoothing contribution is effective 
only within a length scale comparable to the surface diffusion length. 

Fig.5 illustrates the effect of temperature or bond strength and supersaturation on the 
morphology of a growing crystal with the inclusion of both volume and surface diffusion. The 
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geometry and simulation procedures are the same as in fig.4. In fig.5a,b the supersaturation is 
fixed (A|i/kT = 0.69) and 0/kT decreases from 3.9 to 1.6, while in fig.5c,d the temperature 
and bond strength are fixed (<|)/kT = 3.9) and supersaturation increases from 5.0 to 7.0. It can 
be seen that surface roughness can be influenced by either temperature or supersaturation. In 
addition to temperature effects (thermal roughening), surface roughening due to increases in 
supersaturation is evident through the reduction in anisotropy [12]. This leads to growth 
shapes increasingly influenced by volume diffusion. 

Morphological changes under the influence of a single screw dislocation are illustrated 
by fig.6. Besides the dislocation, the geometry and boundary conditions are the same as for 
figs.4 and 5. In order to emphasize the effect of the dislocation the attachment of isolated 
particles (with one solid bond) was suppressed. This approximation is only valid at very low 
temperatures where the attachment probability for the isolated particle is negligible [13]. For 
this figure the supersaturation is fixed (A|i/kT=0.69), and the values of 4>/kT is 3.0, 3.9 and 
5.3, respectively in figs.6a-c. One sees that with a decrease in temperature not only do the 
growth steps become smoother, but also the terrace width increases. Perhaps the feature 
deviating most from the classical theory [21] and results of earlier MC simulations [22,23], is 
this variation of terrace width on the crystal. For this simulation where |a|=b, narrower widths 
are found near the center of the spiral. We explain this as a combination effect of volume and 
surface diffusion. A protrusion is formed at the center of the dislocation during growth. Such a 
protrusion is better supplied by volume diffusion. In addition, due to small overlaps of surface 
diffusion fields farther from the dislocation center, the terrace widths will increase as the step 
velocity increases. Note that in reality, such terrace widening can not be expected to occur in 
vapor systems where |aj >>b 

For fig. 7, a pair of screw dislocations with opposite sign was assumed. As in fig.6, 
the supersaturation was kept at a constant value (A|i/kT=0.69) and <t>/kT was taken to be 3.0, 
3.9 and 4.6, respectively in fig.7a-c. To reveal the competition between 2DN growth at the 
comers and the dislocation spiral growth at the center, the attachment of isolated particles is 
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allowed. Fig.7a shows that at a higher temperature, in spite of a dislocation at the face center, 
growth occurs essentially from the comers of a crystal through 2DN. As the temperature 
decreases or bond strength increases, the dislocation growth at the face center becomes 
increasingly important. At 0/kT = 4.6 (fig.7c) growth occurs only through attachment onto 
steps that originate at the central dislocation pair. Although a supersaturation difference exists 
between comer and face center [12], it is still not high enough to overcome the nucleation 
barrier for 2DN. Instead, the higher supersaturation makes the closed dislocation loop 
unstable, leading to protrusions at comers. Such star shaped dislocation growth has been 
observed experimentally; see, for example, fig.60 in [24]. 
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Figure captions 

Figure 1 Effect of mean free path on growth morphology for 2D square crystal at 

Ap/kT=0.69 and <])/kT=0.69: (a) |a| = 1; (b) |a| = 5; (c) |a|= 10; (d) |a| = 25; (e) 
ty = 50; (f) |a| =100 lattice units. Contours correspond to addition of 1000 
particles each in (a) and 2500 particles in (b)-(f). 

Figure 2 Effect of mean free path on growth morphology for 3D cubic crystal at 
Ap/kT = 0.69 and <j>/kT=2.3: (a)-(b) |a| = 1; (c)-(d) |a| = 5 lattice units. 

Figure 3 Dependence of critical size on mean free path at Ap/kT = 0.69 in both 2D and 
3D for various <J)/kT. 

Figure 4 Effect of surface diffusion on surface morphologies at Ap/kT = 0.69 and 

<j>/kT=1.6: (a) without surface diffusion; (b) with surface diffusion. No volume 

diffusion is considered in both cases. 

Figure 5 Effect of temperature or bond strength and supersaturation on surface 

morphologies at M = lb: (a) A(i/kT = 0.69 and <|>/kT=3.9; (b) Ap/kT = 0.69 and 
<J)/kT=1.6; (c)Ap/kT = 5.0 and <(>/kT=3.9; (d) Ap/kT = 7.0 and 4>/kT=3.9. 

Figure 6 Growth morphologies in the presence of a single screw dislocation with the 

consideration of both surface and volume diffusions at Ap/kT=0.69: (a) <t»/kT= 
3.0; (b) <|>/kT=3.9; (c) <j)/kT=5.3. The growth of isolated particles is suppressed. 

Figure 7 Growth morphologies in the presence of a pair of screw dislocations with 
opposite sign at Ap/kT=0.69. Both surface and volume diffusion are 
included, (a) <t>/kT= 3.0; (b) <j)/kT=3.9; (c) <|>/kT=4.6. 
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We have expanded our earlier Monte Carlo model [Phys. Rev. A38 (1988), 2447; J. Crystal 
Growth 100 (1990), 313] to three dimensions, included re-evaporation after accommodation, 
and growth on dislocation-induced steps. We found again that for a given set of growth 
parameters, the critical size, beyond which a crystal cannot retain its macroscopically facetted 
shape, scales linearly with the mean free path in the vapor. However, the 3-D systems show 
increased shape stability compared to corresponding 2-D cases. Extrapolation of the model 
results to mean free path conditions used in morphological stability experiments leads to order- 
of-magnitude agreement of the predicted critical size with experimental findings. The stability 
region for macroscopically smooth (facetted) surfaces in the parameter space of temperature and 
supersaturation depends on both the surface and bulk diffusion. While surface diffusion is seen 
to smooth the growth morphology on the scale of the surface diffusion length, bulk diffusion is 
always destabilizing. The atomic surface roughness increases with increase in growth 
temperature and supersaturation. That is, the tendency of surface kinetics anisotropies to stabilize 
the growth shape is reduced through thermal and kinetic roughening. It is also found that the 
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solid-on-solid assumption, which can be advantageously used at low temperatures and 
supersaturations, is insufficient to describe the growth dynamics of atomically rough interfaces 
where bulk diffusion governs the process. 

For surfaces with an emerging screw dislocation we find that the spiral growth 
mechanism dominates at low temperatures and supersaturations. The polygonization of a growth 
spiral decreases with increasing temperature or supersaturation. When the mean free path in the 
nutrient is comparable to the lattice constant, the combined effect of bulk and surface diffusion 
reduces the terrace width of a growth spiral in its center region. At elevated temperatures and 
supersaturations, 2-D nucleation controlled growth can dominate in comer and edge regions of a 
facet, while the spiral growth mode prevails in its center. Thus, in addition to confirming the 
experimental observation that the critical size of a growing crystal depends on the prevailing 
growth mechanism, we are able to obtain detailed insight into the processes leading to the loss of 
face and facet stability. 


PACS numbers: 05.50.+q, 61.50.Cj, 61.50.Jr, 64.60.Cn, 68.45.-v. 


I. INTRODUCTION 

Since the first application of the Monte Carlo (MC) method to simulations of crystal 

growth by Chernov and Lewis 1 > 2 its use has become widespread. MC studies were concerned 

with equilibrium and growth morphologies of crystals, 3 * ' 5 * * * * surface roughening transitions and 

growth rate dependence on supersaturation and temperature. 10 ’ 11 For reviews see. 12 ’ 13 These 

studies were primarily focused on interfacial kinetics and have occasionally included surface 

diffusion. 11 A solid-on-solid (SOS) restriction (no overhangs) 14 was often invoked. The 

influence of bulk diffusion, i.e., the transport of growth units through the nutrient to the 

interface, has received less attention until very recently. 15 ’ 18 In the event that surface diffusion 

and interfacial kinetics govern the growth morphology, these simplifications are not severe 
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limitations. However, in reality, bulk transport often plays a decisive role in limiting 
morphological stability. 19 ' 26 

Recently, we have studied the morphological evolution of growing crystals in two 
dimensions by considering both bulk transport and anisotropic interface kinetics. 15 - 16 For the 
bulk transport, we employed the diffusion-limited aggregation model of Witten and Sander. 22 In 
the formulation of the interface kinetics we followed largely Gilmer and Bennema's work. 5 - 11 
Through a systematic variation of the simulation parameters (temperature, bond strength and 
supersaturation), the whole range of growth morphologies from fully facetted to side-branched 
dendritic growth was recovered. Our results show that the diffusion in the bulk nutrient and the 
anisotropy in the interface kinetics act morphologically destabilizing and stabilizing, respectively. 
For a given set of simulation parameters and symmetry of the lattice, there is a critical size 
beyond which a crystal cannot retain its stable, macroscopically facetted growth shape. This 
critical size scales linearly with the mean free path in the vapor. 

In our previous studies 15,16 we ignored the fact that after attachment, the particles still 
have a finite probability to leave the attachment site. This approximation is valid when the system 
is far from equilibrium, i.e., the flux of particles impinging on the crystal surface largely exceeds 
that of the evaporating particles. However, close to equilibrium, the rates of impingement and 
evaporation are comparable. Then a more realistic description of crystal growth must account for 
the non-negligible chance that interfacial particles either move (diffuse) on the surface or 
evaporate back into the nutrient at any time after the initial accommodation on the surface. 

The growth of perfect crystals proceeds typically through two-dimension-nucleation 28 
which is extremely slow at low temperatures and supersaturations. Yet, in reality most crystals 
posses considerable concentrations of defects. Certain defects facilitate growth at low 
temperatures and supersaturations. Frank 29 who found that a screw dislocation can provide a 
inexhaustible source of growth steps. Growth on dislocation-induced steps has since been 
observed experimentally by many workers 30 ' 34 and various theories for this growth mechanism 
have been proposed. 28 - 35 ' 38 However, due to the complexity of the problem, no exact solutions 
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have been obtained. Typical simplifications include the uniform spacing of growth steps, and the 
neglect of nutrient bulk diffusion and surface diffusion. Although there have been some Monte 
Carlo simulations of dislocation-facilitated growth, 39 ’ 40 owing to the computational limitations 
at that time, the transport aspects of the problem were not taken into account. 

In this paper we expand our previous studies ^,16 to three dimension and investigate the 
role of 2-D nucleation growth and dislocation-induced growth in the morphological evolution of 
growing crystals. Instead of the multiple-registration scheme used in the earlier work to reduce 
microscopic noise, 15,16 we employ realistic re-evaporation conditions for interfacial particles. 
The model and simulation procedure are described in Section II. Results on the influence of 
temperature and supersaturation on the surface morphology of perfect crystals are presented in 
Section III.A. Based on a criterion originally proposed by Burton-Cabrera and Frank, 28 the 
thermal and kinetic roughening transitions are quantified. In addition, the scaling between the 
critical stable size of a facetted crystal and the mean free path of the vapor is studied in three 
dimensions. In Section III.B we devote our effort to surfaces with dislocations and pay special 
attention to the shape and spacing of steps of growth spirals. We conclude the paper with results 
on the competition between normal growth at the comers and dislocation-facilitated growth at the 
center of a facet 

II. MODEL AND SIMULATION PROCEDURE 

A. The model 

Two different geometries, planar and spherical, are used to study the effect of bulk 
transport-induced non-uniformity in nutrient distribution on growth morphology. In the planar 
cases (Fig. la), the nutrient is contained between parallel crystal and source planes, which are 
infinite in lateral extent. In order to save computational time, we only consider a portion of the 
whole system and apply periodic boundary conditions in the x- and y-directions. In the other 
cases (Fig. lb) a cubic seed is located at the center of a spherical source. Ideally, the separation 
between the growing crystal and the source should be infinite. But, again, to save computational 
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time we choose a finite separation, yet wide enough that the growth shape remains unbiased.^ 
As a crystal grows, the distance between crystal surface and source is kept constant. 

At molecular length scales, an individual growth unit undergoes generally several basic 
processes before becoming part of the growing crystal. After detachment from the source, a 
growth unit is transported by bulk diffusion towards the growing interface. The actual transport 
kinetics is determined by the interaction of growth units (or their precursors) and other species 
that form the nutrient. Usually, when a growth unit reaches the interface, it is not immediately 
incorporated into the growing crystal. It will adsorb and diffuse on the interface in an attempt to 
find an energetically favorable "final" attachment site, or it will even return to the nutrient before 
it finds such a site. The latter happens, for instance, when the growth unit is misoriented and 
cannot form stable bonds with the crystal. Alternatively, the growth unit may impinge onto a site 
with too few neighbors to prevent it from being dislodged by thermal vibrations before it 
becomes adsorbed. A growth unit may even become dislodged after having arrived at an 
energetically favorable site. The probability that any of these processes occur is determined by 
the local configuration (arrangement of bond-forming neighbors) of the interface sites that the 
unit happens to visit. Though a complete model of such complex scenarios is not practical at this 
point, we have formulated a Monte Carlo model which retains the essential physics of both 
nutrient bulk transport and interface kinetics, including surface diffusion and re-evaporation, 
without resorting to the SOS approximation. 

We assume that the gaseous nutrient phase consists of two components: a growth species 
A, highly diluted in an inert gas B, such that the B concentration is essentially uniform and A-A 
interactions can be ignored. Component B randomizes the motion of A. Convection is ignored. 
The chemical potential is taken to be a linear function of the growth species concentration only. 
This results in a Fickian transport equation and, thus, a random walk may be used to describe 
diffusion in the nutrient phase in the formal® 

U(r,st) = i-X U[r+a, (s-l)x] 


( 1 ) 
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where U(r,sx) is the probability that a walker can be found at location r after s steps (with a time 
interval x) of jump length (mean free path) |aj. The normalization parameter c represents the total 

number of the possible jump sites. 

To approximate the complex processes following the arrival of a growth unit on the 
interface, we make the following assumptions: The impingement rate K + can be obtained, based 
on ideal gas kinetics, from the chemical potential difference Ap. between interfacial vapor and 

average surface site in the form ^ 

K + = Keq exp( Ap / kT) , (2) 

where Keq is the temperature dependent equilibrium value of K + . Note that in reality the overall 
driving force for the diffusion of A towards the crystal and subsequent interfacial attachment is 
the difference in chemical potential between the crystal surface and the source. Since we have 
accounted for bulk diffusion via a random walk process, the chemical potential of the source 
does not appear explicitly in our model. The Ap in Eq.2 is therefor governed by the difference 
between the bulk transport-dependent vapor concentration at the interface and the equilibrium 
concentration at the same temperature. It should not be confused with the chemical potential 
difference (bulk supersaturation, undercooling, etc.) used to control experimental crystal growth. 
This important point is discussed in more detail in 15 - 16 and also below. 

The rates of both the evaporation and surface diffusion processes are sensitive to the local 
configuration of the site from which a unit is to be dislodged. Hence, following Gilmer and 
Bennema, 11 we write the evaporation rate Kj in the site-dependent form 


Kj = v exp(- Ei/kT) 


(3) 
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where v is a lattice vibration factor and, in a nearest neighbor approximation, Ei is simply the 
product of the pair interaction (bond) energy <)) of a unit with a nearest neighbor and n^ the 
number of occupied neighbor sites of site i. 

For surface diffusion we assume that the diffusion rate depends on the occupation 
condition of both the site i that the particle occupies and the potential jump site j. Hence, we 
express the jump rate as 

K i—> j = v s exp( -AEy / kT) f (4) 

where v s is a surface vibration factor, and the activation energy 

r 5i + <(> (ni - n j) for n i > n j 

AE ij= \ (5) 

v 5i for nj < nj 

In reality the term 8j depends on the the specific configuration i-j. For simplicity we have 
assumed that for a given i-j the 5 term is independent of jump direction, hence the notation 5j. As 
a consequence of this assumed direction independence, the Si terms cancel in the following 
formulations of transition probabilities and, thus do not have to calculated. 

In real situations, impingement, surface diffusion and evaporation take place at the same 
time with different magnitude. But in the MC simulation we consider one event at a time. 
Therefore, it is necessary to determine the sequence in which the events are to be considered. To 
this end we define an overall evaporation probability as 

P'=— , (6) 

K + + K* 


where K‘ is the average "evaporation rate 1 
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in 

K ‘=m£ K i' 


i=l 


(7) 


and the summation is over all m interfacial particles. By considering a local equilibrium condition 
for the kink site, on a simple cubic lattice Eq. 6 becomes 

P= E , (8) 

B +exp(A|i/kT) exp(-3<}>/kT) 

with B = — X exp(-i<t> / kT) . 

m i-i 

At equilibrium (Ap/kT = 0), Eq. 8 yields P- * 0.5, i.e. there are approximately equal amounts of 

particles evaporating from and impinging onto the crystal surface. At positive supersaturations 
(Ap./kT > 0), P" is less than 0.5, i.e. impingement exceeds evaporation. Correspondingly, for 
negative supersaturations (Ap/kT < 0) the P" will be larger than 0.5, and more interfacial particles 
will evaporate and the crystal will shrink. These effects will not only be manifested in the growth 
and shrinkage rate, respectively, but will also influence the surface morphology. 

After having decided the sequence of events, the individual probability for the various 
events needs to be determined. This depends on the specific site considered. For impingement, 
the probability, according to Eqs. 2 and 3, is Pi s = K + /(K + + K-). For a simple cubic crystal, 
this can be written as (for a derivation in 2-D see 15 ) 

p s _ exp(Ap/kT) exp[-(3-nj) <>/kT] < ( 9 ) 

l+exp(Ap/kT) exp[-(3-n;) <J)/kT] 

The jump probability from site i to a neighboring unoccupied site j (on the surface or in the 
nutrient), following, 1 5 is 
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p. K i-»j expKnj-nj) tj>/kT)] q 0 ) 

X K i->j X exp[-(n r nj) <j)/kT)] 
j=i j=i 

where c' is the number of unoccupied nearest neighbor sites of site i. Clearly, a larger nj results 
in a higher probability that a molecule will jump to site j on the interface. 

The probability of evaporation for interfacial particles can be expressed as 

e = Kj exp(-n t (fr/kT) i = 1, 2,...m. (1 1) 

i m mB 

Ik; 

i=l 

It can be seen from Eq. 11 that the probability of evaporation decreases exponentially with 
increasing number of solid bonds and, thus, is highly anisotropic. The most probable sites for 
evaporation are those with the least solid neighbors, such as admolecules with a single solid 
bond. As the temperature and, thus, the surface roughness is increased, the anisotropy of the 
evaporation is decreased. 

B. Simulation procedure 

A summary of the simulation steps is presented in Fig.2. First the initial conditions are 
chosen: for growth onto perfect crystals a smooth surface is set at a certain distance (typically 50 
lattice units b, for |a| =b) from the source. In the planar case the surface is chosen as a square 
lattice measuring 60b x 60b, with periodic boundary conditions in x and y (lateral) directions. 
For growth onto dislocated crystals the lateral size is doubled in order to provide space for the 
evolution of several turns of a spiral. An initial dislocation is introduced by a vertical slip of parts 
of the lattice similar to. 40 Note also that, as the crystal grows, the source is moved outwards 
such that the separation remains unchanged, as discussed in. 15 

Following the setting of initial conditions and choice of specific input parameter values, 
P- is calculated from <}>, T, Ap. and the current crystal morphology. Comparison of a random 
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number R with P' determines the event to be considered. If R > P', a growth unit is released 
from the source and bulk diffusion is simulated through an isotropic sequence of random jumps 
of equal length (mean free path) |a|. Only when the particle has come within a distance |a| from the 
crystal surface is a check implemented that determines whether the interface has been reached. 
The following steps, including the determination of average surface diffusion lengths, are similar 
to our previous simulations except that we do not impose the multiple registration used in. 1 ^’ 

A growth unit is tracked until it has either escaped from the system or is stuck onto the crystal 
surface. 

If R < P", an evaporation event is selected. Then, all Pi e (i=l,...m) are calculated from 
the recorded positions and neighboring configurations. Another random number R' is generated. 
If by chance R' falls into i 1 ^ interval, the i 1 ^ 1 interfacial particle is chosen to move. Then, 
governed by Eqs. (9), (10) and (11), the particle is kept moving until it has attached to some 
interfacial site or escaped to the source. If it sticks on site j, the local Pj e s, around the new site j 
as well as the earlier site i where the particle originated, are recalculated. By using a new random 
number, a new evaporation event is selected. The evaporation route is pursued until some particle 
has escaped to the source. When this happens the P* is recalculated and a new event is selected. 

The simulation is continued until a desired size or layer thickness is reached. In a planar 
case, the simulation is continued until 14,400 or 10,000 particles, respectively, have been added 
to a perfect or dislocated surface. In a spherical case, the simulation is stopped after a crystal has 
become unstable (the criterion for the instability will be discussed later). Typical CPU times were 
around 3,000-10,000s on a CRAY X-MP/24 and 5 - 50 hours on an ARDENT computer. 

Our model has been designed for computational efficiency. As outlined above (see Eqs. 
6-8), the event to be considered is decided upon before the execution based on physical criteria 
rather than random choice only. In this way, any "unnecessary" operation is avoided. For 
example, evaporation will be executed whenever a random number R is less than P*, although 
the selection of the actual execution site is quite random. In the earlier models 11(40 the selection 
of a particle and calculation of the evaporation probability do not necessarily lead to an actual 


execution of an evaporation process. Many interfacial particles may have to be selected, before a 
real evaporation event takes place. This is particularly time-consuming under conditions of low 
temperature and supersaturation, in which the interface is rather smooth and the probability of 
evaporation is very low. 

Note that, unless stated otherwise, the simulations were carried out with the magnitude of 
the mean free path in the vapor equal one lattice constant, lal = b. 

III. RESULTS AND DISCUSSION 

A. Perfect surfaces 

1. Effects of surface and bulk diffusion at various temperatures 

As references for the more complex cases to be treated later and for comparison with 
earlier work, 6 ' 9 we have first examined cases in which both bulk and surface diffusion were 
neglected. Particles arriving at random locations of the crystal surface either stick or are 
discarded. Figure 3 shows the effect of growth temperature or bond strength (i.e. 4>/kT) on the 
morphology of a planar crystal surfaces at a constant low value of Ap. /kT = 0.35. At low 
temperature, as can be seen from Fig.3a, the growing crystal is atomically smooth with 
occasional 2-D clusters on the surface, i.e. growth occurs by 2D nucleation only. 28 As the 
temperature increases (Fig. 3b) the surface becomes rougher and more 2D clusters with smaller 
size are formed on the surface and even on top of large clusters. On further increase in 
temperature the whole surface becomes atomically rough (Fig.3d), and large size clusters can 
hardly be discerned. In this situation the crystal grows more readily due to the disappearance of 
the energy barrier for 2D nucleation. This thermal roughening effect has been extensively studied 
before under equilibrium conditions. 6 * 9 -* 2 ’^ The earlier studies basically employed an algorithm 
for random pair exchange controlled by a Boltzmann factor exp(-AE/kT). At equilibrium, the 
validity of the algorithm is guaranteed by the thermodynamic principle of path independence. For 
a non-equilibrium situation the random pair exchange algorithm is no longer valid. The final 
result depends not only on the potential energy between the initial and final sites but also on the 


actual path of the particles. In contrast, our model, as discussed in section II, is more general and 
can be used in both equilibrium and non-equilibrium situations. 

To further demonstrate the effect of temperature on the morphology of a growing crystal 
we have plotted in Fig.4 the surface area, normalized by the area of a perfectly smooth surface, 
versus number of attached particles (in units of completed layers). The four curves correspond to 
cases a-d in Fig. 3. As can be seen from Fig.4, at the lowest temperature the surface area 
oscillates periodically with minima on layer completion, alternating with maxima at about half- 
filled layers. Thus growth proceeds essentially layer by layer. When, following the formation of 
a 2-D nucleus, a layer spreads, the surface area increases until a half-filled layer state is reached 
and then it decreases until a layer is completed. However, when the temperature is increased and, 
thus, the surface becomes rougher, the distinction between a half-filled layer and completed 
layers is blurred, i.e. new layers form before completion of the earlier layers. With increasing 
roughness, the periodicity disappears and the surface becomes delocalized. 

Figure 5 illustrates the consequences of surface diffusion. Bulk diffusion is not 
considered in these examples with 4>/kT= 1.6, A[i/kT=0.69 . It can be seen that with surface 
diffusion (Fig. 5b) the crystal surface is much smoother than without (Fig.5a). The reason is that 
surface diffusion can provide an additional way for interfacial particles to relax to some 
energetically more favorable (low energy) sites and, hence, reduces the amplitude of the surface 
roughness. Of course, as we have shown before, 16 such a stabilizing contribution is effective 
only at a length scale comparable to the surface diffusion length, which, in turn, is sensitive to 
the surface roughness. 16 Generally speaking, the smoother a surface is, the farther an interfacial 
particle can diffuse during its lifetime on the surface. Hence, as shown by the simulation results 
of Fig.6, the surface diffusion length decreases with increasing temperature. In addition to the 
temperature effect, Fig.6 also reveals a strong dependence of surface diffusion length on 
supersaturation. This was not taken into account in earlier simulations, 11 in which the surface 
diffusion length was considered as an externally adjustable parameter, independent of 



temperature and supersaturation. By contrast, in our model surface diffusion is an integral part of 
the attachment and evaporation processes. 

Unlike surface diffusion, bulk diffusion acts destabilizing on the interface morphology. 
15-17 This is shown in Fig.7 for a case with surface diffusion, A|i/kT =3.0 and <J)/kT =2.3. For 
the three subfigures the mean free path in the nutrient was decreased from "infinite" (Fig.7a), to 
ten (Fig.7b) and one crystal lattice constant (Fig.7c), respectively. The infinite mean free path 
corresponds to direct jumps of growth units from the source onto the crystal surface, as in the 
simulations leading to Figs. 3 and 5. A mean free path of one lattice unit implies an extremely 
dense nutrient similar to a liquid. It is evident from this figure that bulk diffusion, particularly in 
the case of a short mean free path (Fig.7c), destabilizes the surface morphology. With increasing 
surface roughness not only is the stabilizing effect of surface diffusion weak, but simultaneously 
the stabilizing anisotropy in surface kinetics is decreased. Once a protrusion forms by chance, it 
will be amplified by the bulk diffusion. At intermediate values of temperature and 
supersaturation, there exists a competition between surface kinetics and bulk diffusion. Only if 
the surface kinetics prevails, can crystal surfaces remain stable, i.e. facetted. 16 More pronounced 
destabilization from bulk diffusion at higher temperatures will be shown later. 

2. Supersaturation effects 

Increases in supersaturation can also affect the crystal morphology through kinetic 
roughening. 41 " 44 This is also born out by our 2-D simulations. 15 ' 16 The influence of 
supersaturation on a 3D growth morphology is shown in Fig.8, in which the temperature and 
bond strength are kept constant (<t>/kT = 3.9). Both surface and bulk diffusion are included. One 
sees that with increasing supersaturation the size of the 2-D clusters decreases, the surface 
roughens, eventually leading to many depressions and protrusions on the surface, and even 
vacancies in the crystal. As formulated in section II.A, an increase in supersaturation causes an 
increase in impingement rate (Eq.2) and, hence, a relative decrease in evaporation probability 
(Eq.6). When the supersaturation reaches some critical value, the surface kinetics become 
relatively less important and the growth process is controlled by bulk diffusion. The ensuing loss 
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of facet stability is further illustrated in Fig.9 by profiles of the interfacial layer occupation 
numbers. The four profiles correspond to the cases depicted in Fig. 8. One clearly sees that the 
interfacial width grows with increasing supersaturation. 

The above results reveal significant concentrations of overhangs and vacancies in the 
crystal at higher supersaturations and temperatures. Hence, under these conditions, the SOS (no 
overhang) assumption is unrealistic for the description of crystal growth morphologies. Thus, 
the SOS assumption is only a good approximation at low temperatures and supersaturations 
where the surface morphology is controlled by surface kinetics rather than bulk diffusion. Note, 
however, that the supersaturation value required for significant kinetic roughening depends on 
both bond strength and temperature. A crystal surface can have the same roughness at different 
combinations of supersaturation and temperature. 

In order to illustrate the dependence of the anisotropy in interface kinetics parameters on 
growth conditions, we have plotted in Fig. 10 the sticking probabilities at sites with various 
numbers of nearest neighbors as a function of temperature/bond strength and supersaturation 
(Eq. 9). For clarity only surfaces for Pi s , P 3 S (kink site) and Ps s are plotted in this figure. Note 
the temperature independence of P 3 S that results from the assumption of local equilibrium at the 
kink site underlying Eq.9. Figure 10 well illustrates the strong dependence of the anisotropy in 
sticking probability (separation between Pi s and Ps s ) on temperature and super saturation. At low 
values of these parameters, Pi s and Ps s are nearly constant and close to 0 and 1 respectively. 
This reflects the difficulty of the attachment of isolated particles and the ease with which holes 
(i.e. sites with more than 3 solid nearest neighbors) are filled. One can also see from this figure 
that the anisotropy in sticking probability can be more effectively reduced by increasing the 
temperature than by increasing the supersaturation. There is a rapid decrease in the anisotropy of 
the sticking probability around <J)/kT = 1.0 - 2.0 at low supersaturation. The value of the 
anisotropy in sticking probability is a fundamental quantity for controlling the morphology of a 
crystal surface. Only when the anisotropy is significant can a crystal retain a stable, facetted 


form. 
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3. Surface roughness 

To further quantify the surface morphology features we use a criterion originally 
introduced by Burton-Cabrera and Frank (BCF) 28 that defines the surface roughness R s in terms 
of the surface energy E (i.e. the number of broken bonds times <f>) at temperature T and the 

surface energy Eq for a perfectly smooth surface at zero temperature, in the form 



This R s is a merely geometrical measure that depends only on the number of broken lateral 
bonds per unit area. The variation of R s with temperature obtained from a simulation of an 
equilibrium situation (A|i/kT = 0) is plotted in Fig. 11. One can see that the surface roughness 
monotonically increases with temperature, with a slope that is highly temperature dependent. 
There is also a unique inflection point (location of maximum slope), which, as better indicated by 
the plot of the derivative dRs/dT(T) in the same figure, lies at kT/<J> = 0.62. Since R s is 

proportional to the surface energy, dRs/dT should be proportional to the heat capacity or specific 
heat of the surface. A singularity in heat capacity is characteristic for a phase transition. But in 
MC simulation, the singularity has been degraded owing to the effect of finite system size. 45 ’ 46 
Swendsen 9 first used a maximum in heat capacity in a MC simulation to define the roughening 
transition temperature. He obtained for a simple cubic lattice Tr = 0.575 <}>/k. Leamy and Gilmer 

7 found, by using a different criterion and employing the SOS restriction, a roughening transition 
temperature Tr => 0.64 <>/k, which is surprisingly close to our value. This indicates again that the 

SOS approximation is appropriate for the description of surface morphology under equilibrium 
conditions. However, as has been shown above, the SOS assumption is no longer valid when 
the growth becomes bulk diffusion controlled, i.e. the surface has become rough. 

In Fig. 12 we have delineated the temperature- supersaturation combinations for surface 
roughening transitions, as obtained from our MC results and Eq.12 with the inflection point 


criterion, for systems with various combinations of surface and bulk diffusion conditions. As 
reflected by the extent of the respective stability regions (smooth vs. rough), surfaces on which 
surface diffusion is significant can retain atomic smoothness up to higher temperatures and/or 
supersaturations than without surface diffusion. Bulk diffusion, on the other hand, reduces the 
smooth (stable) region in the <>/kT vs. Afi/kT plane, as we have seen more qualitatively before. 

4. Re-evaporation effects 

All of the above results were obtained with the new algorithm that accounts for the 
possibility of re-evaporation of interfacial particles. To gauge the effect of re-evaporation with 
respect to the earlier noise reduction algorithm, 15,16 Pig. 13 presents some results obtained 
without re-evaporation and either without or with the multiple-registration scheme used in our 
earlier work. Otherwise the growth conditions in Fig. 13 a-b are the same as in Fig. 8a ( A|i/kT = 
0.69) while the conditions in Fig.l3c-d are the same as in Fig.8c (Ap/kT = 5.0), with <|)/kT = 3.9 
for all. The multiple-registration noise reduction scheme 15 * 16 used for Figs. 13b and 13d requires 
10 registrations of incoming particles at a site i before that site is considered occupied. 
Comparison of the corresponding cases in Figs. 8 and 13 shows that at low supersaturations the 
noise reduction scheme is somewhat more effective in smoothing the surface than re-evaporation 
At high supersaturation, however, both schemes give nearly the same smoothness (compare 
Figs. 8c and 13d). Similar conclusions can be drawn for the smoothing effect of the two schemes 
upon variations of 0/kT. Hence we conclude that at situations far from equilibrium or at high 

temperatures, the noise reduction scheme can give qualitatively similar results as the new 
evaporation scheme, for which, however, the physical implications are much clearer. 

5. Spherical source 

The above results were based on the uniform nutrient flux conditions of the planar 
geometry of Fig. la . But in reality nutrient flux conditions are most often non-uniform. To 
further explore the effect of nonuniform supply beyond the 2-D results obtained in 15 > 16 with a 
circular source, we have performed MC runs with growth onto a cubic surface inside a spherical 
source; see fig lb. Figure 14 presents surface morphologies obtained in this geometry with <t>/kT 
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= 2.3 and A|i/kT = 0.69. The mean free path or jump length is one lattice constant, i.e. |a| = b in 
Fig.l4a-b, and five lattice constants (]a| = 5b) in Fig.l4c-d. Similar to the 2-D findings in, 15,16 
one sees that loss of stability of a facet is associated with the formation of a depression in its 
center and, correspondingly, preferred growth at comers and edges. As discussed before 
this occurs because the anisotropy in surface kinetics can compensate for the nonuniformity in 
nutrient supply only up to a certain critical size of the crystal. This nonuniformity decreases with 
increasing mean free path of the diffusing nutrient particles. Hence, as is shown by Fig. 14, the 
critical crystal size increases with increasing mean free path. To quantify the critical size we have 
plotted in Fig. 15 the normalized total surface area of the growing cubic crystal versus the size of 
the crystal . The surface area is normalized by the total surface area (six faces) 6M 2 ^ 2 of a 
perfectly smooth cubic crystal consisting of M growth units. As long as the surface area follows 
the 6M 2 ^behavior, a crystal is considered morphologically stable. When a crystal begins to lose 
its stable facetted form its normalized total surface area increases. Somewhat arbitrarily, we 
define the critical size as the one at which the surface area exceeds that of a smooth surface by 
15%. Thus we find that the crystals depicted in Fig. 14 have critical sizes of about 27b and 52b, 
respectively , for |a| = b and |a| = 5b. 

In Fig. 16. we have summarized our 3-D results for the dependence of the critical size on 
mean free path at fixed Ap./kT and various 4>/kT 's, together with the earlier 2-D results. 16 One 
can see that the critical size increases with decreasing temperature. At fixed temperature and 
supersaturation, the critical size scales linearly with the mean free path in the parameter range 
considered. The 3-D simulations yield larger slopes than the corresponding 2-D cases, i.e. the 3- 
D cases are morphologically more stable. This results from the fact that, on average, there are 
more solid neighbors associated with interfacial particles in three dimensions and, thus, the 
(stabilizing) anisotropy in interface kinetics is more pronounced. Although computational time 
limitations do not allow for the direct simulation of the critical size at the mean free path used in 
morphological stability experiments, 22 - 26 Fig. 16 encourages us to scale linearly to these 
conditions. This leads to an order-of-magnitude agreement between experimental and modelling 


results for the critical size, with the 3-D results approaching the experimental findings closer than 
the 2-D results. 16 

B. Surfaces with dislocations 

1. Temperature effects 

The development of a growth spiral from an initially straight step that results from a 
screw dislocation with Burgers vector of one lattice constant normal to the (001) face is 
illustrated in Fig. 17. The bond strength/temperature (i.e. <J>/kT) dependence of the surface 
morphology with a single screw dislocation is shown in Fig. 18. In this simulation, both surface 
and bulk diffusion were ignored.The supersaturation was kept constant (Ap/kT=0.69) while the 
<J)/kT was decreased successively from 5.3, 4.6, 3.0 to 2.0. One sees that at low temperature or 
high bond strength (Fig.l8.a) the steps of the resulting growth spiral are quite smooth and the 
shape is highly polygonized, leaving essentially only low-index steps exposed. As the 
temperature is increased (Fig.l8b-d) not only do the steps roughen but also the shape of the 
spiral becomes rounder and eventually indiscernible as 2-D nucleation becomes pronounced. 
These morphological changes are a result of the increase in surface roughness with growth 
temperature. At low temperatures (Fig. 18a), growth occurs only via attachment to the spiral 
steps, since the roughness of the remainder of the facet is too low to result in significant sticking 
probabilities. Furthermore, due to the low roughness of the steps, the step attachment kinetics is 
highly anisotropic, leading to rather straight step shapes. As temperature is increased, the steps 
and remainder of the face roughen. This results in a reduction of the anisotropy in step 
attachment kinects and, thus, rounding of the steps. Simultaneously, with increasing face 
roughness, the energy barrier for 2-D nucleation is reduced. On further temperature increase 2-D 
nucleation-assisted growth becomes increasingly important, until it dominates at 4>/kT = 2.0. 

These results confirm that it is reasonable to neglect 2-D nucleation in simulations of 
dislocation-assisted growth morphologies at low temperature and supersaturation, where the 
sticking probability Pi s (see Fig. 10) is very small. This allows for a drastic reduction in 



computational times required. 39 Hence, throughout the following simulations of dislocation- 
assisted growth, the attachment of isolated particles, i.e. with one bond to the surface, is 
suppressed. 

Figure 19 was obtained by considering a pair of screw dislocation with opposite sign, 

i.e. a Frank-Reed step source. 47 In order to save computer time and memory, the separation 
between the centers of the dislocations was chosen to be only 7b. In reality, this separation may 
be much larger. When the two growth spirals turn in opposite direction and meet, closed loops 
are formed periodically. During this sequence of simulations the supersaturation was kept 
constant while <j)/kT was decreased from 5.3 to 4.6 and 3.9, respectively, and surface and bulk 
diffusion were ignored. Again, as the temperature increases the steps become rougher and the 
closed loops become rounder and spaced more closely, analogous to the temperature dependent 
behavior displayed in Fig. 18. 

2. Supersaturation effects 

The effect of supersaturation on the growth morphology of a face with a single screw 
dislocation is shown in Fig.20, in which the <J)/kT is kept at 5.3 while the A|i/kT is increased 
successively from 0.69 to 2.0 and 3.0, again ignoring surface and bulk diffusion. As can be 
seen, increases in supersaturation also make steps rougher and less polygonized, similar to the 
effect of 4>/kT in Figs. 18 and 19. This morphological change is, of course, due to kinetic 
roughening of the steps, in contrast to the thermal roughening occurring as the temperature is 
increased. In addition, the terrace widths between adjacent spiral arms decrease with increasing 
supersaturation, consistent with the prediction of classical theories.^ 3 ^' 3 ^ 

3. Effects of surface and bulk diffusion 

As has been discussed in section III.A, surface diffusion can greatly smooth the 
morphology of perfect surfaces. The effect of surface diffusion in the presence of a single 
dislocation is demonstrated in Fig.21. For this sequence, nutrient bulk diffusion was ignored 
and the growth conditions (i.e. <J)/kT and A|i/kT) were the same as those for Fig.20. We find, in 
comparison to Fig.20, not only smoother steps, but also more pronounced polygonization. This 
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is because with surface diffusion, more growth units can reach rounded comers, where, due to 
the higher kink density, the attachment probability is higher. The spirals also become more 
polygonized due to the resulting increase in anisotropy of step attachment kinetics This is at 
variance with the conclusion of Sunagawa and Bennema, 3 ^ who expected that surface diffusion 
will suppress the formation of close-packed or periodic-bond-chain oriented steps. 

Growth morphologies of a dislocated surface in the presence of both surface and bulk 
diffusion are depicted in Fig. 22. The growth conditions are the same as in Figs. 20 and 21, i.e. 
<j)/kT is kept at 5.3 and the supersaturation is increased. It is evident from this figure that bulk 
diffusion is also destabilizing for steps, as revealed, for instance, by the lateral depressions at 
high supersaturation (Fig.22c). But probably the most interesting new feature revealed by this 
simulation is the variation of the terrace width, with narrower terraces found near the center of 
the spiral, particularly at higher supersaturation (Fig.22c). This is at variance with earlier MC 
simulations without bulk diffusion 3941 and classical theories. 8 . 35 ' 38 Realizing that the above 
simulation is based on a mean free path length (in the bulk) equal to a lattice constant, the 
decrease in terrace width toward the center of the spiral can be understood as a combined effect 
of bulk and surface diffusion. Since the growth hillock protrudes into the nutrient, steps and 
terraces close to its center of the spiral are somewhat better supplied with growth units by bulk 
diffusion. This leads to tighter winding of the spiral.^ 8 On the other hand, the surface diffusion 
fields of terraces at the periphery of the spiral overlap less than those near the center. Overall, 
this leads to higher spreading velocities of the outer turns of the spiral and, thus, an increase in 
terrace width with distance from the spiral's center. In real vapor systems, where |a| »b, such 

behavior should not be expected. It is more likely to occur in growth from condensed phases. 

4. Spherical source 

All of the above results for dislocation-assisted growth were, obtained for the planar 
geometry of Fig. la. To demonstrate the competing effect of 2-D nucleation growth and 
dislocation growth in a non-uniform concentration field, we have carried out a simulation for the 
3-D geometry of Fig.lb with the same growth conditions as for Fig.l4b. Among the six faces of 
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the cubic crystal, five are assumed to be perfect and only one has a pair of dislocations with 
opposite sign. Surface diffusion is taken into account. Figure 23 shows the crystal at a stage 
after exceeding its critical stable size (see discussion in section m.A). The subfigures show 
different faces of the same crystal, with (a) displaying three of the five faces without dislocation 
and (b) the dislocated face on top. The dislocation pair is seen to enhance the face stability. All 
five originally perfect faces have developed central depressions similar to that found in Fig.23a. 
This finding can be understood in terms of the different supersaturation dependence of the 
growth rate in 2-D nucleation and dislocation-assisted growth. Bulk diffusion results in a lower 
supersaturation in the face center than at the edges and comers. Yet, the face center region, due 
to the presence of dislocation-induced surface steps, offers a higher probability for the 
attachment of growth units than the edge regions in which steps are generated only by 2-D 
nucleation. This increased attachment probability in the center of the face can compensate for the 
leaner supply of growth units into this region, see also.^.N* This competition of growth steps 
originating at dislocations in center regions of facets, with growth steps from 2-D nucleation near 
comers, has been experimentally observed by several workers; for references see.^8 Also, the 
better utilization of growth units (through the anisotropy of kinetics coefficients) in regions that 
are less readily supplied by the bulk diffusion field is the key point of Chernov's anisotropic 
stability theory for facet growth. 21 

To further illustrate the competition between 2-D nucleation growth and dislocation 
growth in a non-uniform nutrient field, we have carried out a simulation at a much lower 
temperature (<J>/kT = 5.3) than that used for Fig.23. At the lower temperature the surface 
roughness and, thus, the sticking probability is considerably reduced. Hence, in order to not 
exceed prudent computational times, we have modelled only the 3-D evolution on one face, with 
special periodic conditions on the lateral boundaries of the pyramid outlined in Fig. lb. Note that 
the outer boundary condition in this configuration still corresponds to a spherical source as in 
Fig.23. In this simulation, <j)/kT is kept at 5.3, and the supersaturation is increased from 0.6 to 
2.0 and 3.0, respectively. The simulations include surface diffusion and unsuppressed 
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attachment of single particles. Figure 24 shows the results after attachment of 20,000 particles in 
all three cases. One sees that at low supersaturation (Fig.24a) growth occurs only through 
attachment onto steps that originate at the central dislocation pair. Under these growth condition, 
the supersaturation increase at the comers 16 is not high enough to overcome the nucleation 
barrier for 2-D nucleation, which, due to the low surface roughness, is high. The higher 
supersaturation at the comer, though, destabilizes the dislocation-induced growth step loops, 
leading to lateral protrusions towards the comers, similar to the growth patterns obtained in the 
2-D situation. 16 ’ 16 Such star shaped dislocation growth has been observed experimentally; see, 
for example, Fig.60 in. 34 As the supersaturation is increased, the controlling effect of the central 
dislocations decreases and 2-D nucleation becomes significant at the comers (Fig.24b). At an 
even higher supersaturation (Fig. 24c), growth is essentially dominated by 2-D nucleation at the 
comers, in spite of the dislocations at the face center. 

IV. CLOSING REMARKS 

The above simulations, as well as most earlier efforts to model surface morphologies, are 
based on the supersaturation/temperature and bond strength/temperature parameters Ap/kT and 
respectively. These are highly idealized scaling parameters, which require utmost caution 
in attempts to quantitatively compare the model predictions with actual experiments. 

Specifically, it must be reemphasized that, as discussed in section II. A in connection with 
Eq.2, the chemical potential difference or supersaturation used in these models is not equal to the 
bulk nutrient supersaturation typically determined by the experimentalist. This chemical potential 
difference is solely that part of the overall difference that drives the attachment of growth units 
once they have been transported to the interfacial region, say to within a mean free path of the 
growing surface. Though possible in principle, no unambiguous measurements of this interfacial 
supersaturation have become available as yet. However, all crystal growth theories, as well as 
measurements of interfacial undercoolings in facetted and non-facetted regions of an interface that 
grows from a melt 49,50 indicate that this interfacial chemical potential difference necessary to 


23 


drive a certain attachment (growth) rate depends on the local interface morphology (kink and step 
density, etc.) and, thus, on the locally governing growth mechanism. This fact is ignored by all 
modelling at this point by fixing the interfacial A p irrespective of the "underlying" surface 

morphology. 

With respect to correlations of the bond or pair interaction energies, used in these models, 
with values of actual systems, similar caution is required. In addition to the fact that the simple, 
highly symmetric bond picture is hardly an accurate representation for actual atomic (quantum 
mechanical) interactions, one must realize that the magnitudes for bond strengths are traditionally 
derived from averaged bulk properties, rather than the (relaxed) surface states that govern the 
attachment kinetics. Some, small improvement of this coarse description, based on a "variable 
bond model", has been made by our group earlier. 51 

The above material shows that current kinetic crystal growth modelling can provide 
considerable physical insight into the effect of various growth parameters. Yet, it is also clear 
that much more experimental and theoretical work is required before quantitative fidelity can be 
expected. 
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Figure captions 

Fig. 1. Schematics of simulation geometries, (a) Planar case: lateral dimensions 60 x 60 

lattice constants, with periodic boundary condition in the x- and y-directions. Space 
between crystal surface and source contains nutrient through which growth units 
diffuse to the crystal surface, (b) Spherical case: growing crystal in center of spherical 
source. Pyramidal section is used with periodic boundary conditions in azimuthal 
direction for simulation of 3-D growth morphology at low temperatures; see text. 

Fig. 2. Flowchart of simulation steps. 

Fig. 3. Effect of temperature and bond strength on crystal surface morphologies at fixed 
supersaturation, A(i / kT = 0.35. Both surface and bulk diffusion are excluded. 

(a) <J>/kT=3.9, (b) <j>/kT=2.0, (c) <J>/kT=1.6, (d) <J>/kT=1.4. 

Pig. 4. Normalized surface area as function of number of attached particles (in units of full 
layer, i.e., 3600 particles). The normalization is based on the area of a perfectly 
smooth surface. The growth conditions are the same as in Fig.3. 

Fig. 5. Effect of surface diffusion on surface morphologies at <J> / kT = 1.6 and A|i / kT= 

0.69 without consideration of bulk diffusion, (a) without surface diffusion, (b) with 
surface diffusion. 

Fig. 6. Dependence of average surface diffusion length on bond strength/temperature at two 
different supersaturations, Ap / kT= 0.69 and 3.0. 

Fig. 7. Effect of the mean free path length in bulk diffusion on surface morphologies at 

(j) / kT = 2.3 and Ap / kT = 3.0. Surface diffusion is included, (a) |a|=«> (b) |a| =10b, 

(c) |a|=lb. 

Fig. 8. Effect of supersaturation on surface morphologies at <J> / kT = 3.9 and a mean free 
path of one lattice unit, (a) A|i/kT=0.69, (b) Ap/kT=3.0, (c) Ap/kT=5.0, (d) 
Ap/kT=7.0. 

Fig. 9. Coverage of interfacial layers in the surface morphologies of Fig.8. 
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Fig. 10. Sticking probability Pj s (Eq. 9) as a function of <]>/kT and Ap/kT. For clarity, only 
Pl s , P 3 S and P 5 S are shown. 

Fig. 1 1 . Surface roughness and its derivative with respect to temperature as function of <|)/kT at 
equilibrium (Ap/kT=0). The numbered arrows on the abscissa indicate roughening 

temperatures obtained earlier by . ' ° 

Fig. 12. Stability diagram for smooth (stable) and rough (unstable) growth morphologies as 
a function of <>/kT and Ap/kT. SD: with surface diffusion; BD: with bulk diffusion. 

Fig. 13. Effect of noise reduction (multiple registration)^ on surface morphologies at 

<j>/kT = 3.9. (a) Ap/kT = 0.69, without noise reduction; (b) Ap/kT = 0.69 with 
noise reduction; (c) Ap/kT = 5.0, without noise reduction; (d) Ap/kT = 5.0, with 
noise reduction. 

Fig. 14. Growth morphology of three dimensional cubic crystal in a spherical source at 4>/kT 
= 2.3 and Ap/kT = 0.69. Noise reduction through multiple-registration. (a,b) |a| 

=lb, (c, d) |a| =5b.(a) On addition of 6340 particles to nucleus with 21x21x21; (b) 
Further addition of 3 1056 particles to (a); (c) On addition of 28828 particles to 
nucleus with 47x47x47; (d) Further addition of 52542 particles to (c) 

Fig. 15. Normalized total surface area as function of crystal size for a growing cubic crystal 
The growth conditions are the same as in Fig. 13a and b. 

Fig. 16. Dependence of critical size on mean free path at Ap/kT = 0.69 and various <|>/kT for 

both 2D 16 and 3D results. 

Fig. 17. Time evolution of growth step originating at single screw dislocation .<}>/kT = 5.0 and 
Ap/kT=0.69. (a) initial straight step, (b) after attachment of 1000 particles, (c) after 
attachment of 5000 particles. 

Fig. 18. Bond strength/temperature dependence of surface morphology of face with a single 
dislocation. Ap/kT=0.69. Both surface and bulk diffusion are ignored. 

(a) <J)/kT = 5.3, (b) <{>/kT = 4.6, (c) <>/kT = 3.0, (d) <j>/kT = 2.0. 
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Fig. 19. Bond strength/temperature dependence of surface morphology of face with a pair of 
dislocations with opposite sign. A(j/kT=0.69. Attachment of isolated particles 
suppressed, both surface and bulk diffusion are ignored, (a) <)>/kT = 5.3, 

(b) <j)/kT = 4.6, (c) <t>/kT = 3.9. 

Fig. 20. Supersaturation dependence of surface morphology with a single screw dislocation. 
<j)/kT = 5.3. Both surface and bulk diffusion are ignored, (a) A|i/kT=0.69, 

(b) A|i/kT=2.0, (c) Ap/kT=3.0. 

Fig. 21. Effect of surface diffusion on surface morphology with a single screw dislocation. 

Same growth conditions as in Fig. 20. Bulk diffusion is neglected. 

Fig. 22. Effect of surface and bulk diffusion on surface morphology with a single screw 
dislocation with same growth conditions as in Fig.20. 

Fig. 23. Growth morphology of a cubic crystal in a spherical source. The growth conditions are 
the same as in Fig. 14b. Different faces on the same crystal, (a) all faces without 
defects, (b) top face with pair of dislocations of opposite sign. 

Fig. 24. Growth morphology of a dislocated face in a spherical source at <J)/kT = 5.3 with bulk 

and surface diffusion, and unsuppressed formation of 2D nuclei. 

(a) A|x/kT=0.69, (b) A^i/kT=2.0, (c) Ap/kT=3.0. 
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ABSTRACT 

A Monte Carlo model is used to simulate the morphological evolution of crystals growing 
from an incongruent vapor phase. The model combines nutrient transport, based on a modified 
diffusion- limited aggregation process, with anisotropic surface kinetics and surface diffusion. 
Through a systematic variation of the simulation parameters (temperature, bond strength and 
supersaturation), the whole range of growth morphologies from fully facetted to side-branched 
dendritic growth is recovered. 

The diffusion in the bulk nutrient and the anisotropy in the interface kinetics are seen to be 
morphologically destabilizing and stabilizing, respectively. It is found that for given bond strength 
and symmetry of the lattice, there is a critical size for stable, macroscopically facetted growth. This 
critical size scales linearly with the mean free path in the vapor. Since both thermal and kinetic 
roughening reduce the kinetic anisotropy, the critical size decreases as either temperature or 
supersaturation are increased. Surface diffusion is seen to stabilize facetted growth on the shorter 
scale of the mean surface diffusion length. 

In simulations with a uniform drift superimposed on the random walk nutrient transport, 
crystal faces oriented towards the drift exhibit enhanced morphological stability in comparison to the 
purely diffusive situation. Rotational drifts with periodic reversal of direction are found to be 
morphologically stabilizing for all faces of the crystal. 

INTRODUCTION 

The characterization of the conditions under which a growing crystal is capable of preserving 
its shape, i.e. is morphologically stable, is both scientifically challenging and technologically 
important. Morphological stability is necessary for the growth of homogeneous single crystals that 
are needed for numerous device applications. The growth morphology and morphological stability of 



a crystal, however, is an extremely complex problem, involving, in general, an interaction between 
nutrient transport and interface kinetics [1-3]. 

Two basic approaches to the morphological description of crystal growth have been adopted. 
The first is the macroscopic approach which involves the solution of a continuum transport equation 
coupled with moving boundary conditions [4], The second is the simulation of microscopic 
processes by tracking the individual growth units [5]. 

The continuum approach has lead to significant insight into the morphological evolution of 
essentially isotropically responding, i. e. nonfacetted interfaces, such as prevail in many melt growth 
systems The theory of morphological stability of non-facetted growth is well established. Stability 
conditions are known for a variety of simple geometries [4,6-9], including systems with weakly 
nonlinear [9-11] and highly nonlinear isotropic response [12-14]. In addition, the coupling between 
morphological and hydrodynamic instabilities has been examined [15-20]. Continuum models have 
also been developed for the morphological stability of nonfacetted interfaces with anisotropies in 
surface tension and growth kinetics [e.g. 21-24], One of the most important results of these 
treatments is the insight that anisotropies can stabilize otherwise unstable closed growth forms up to 
a critical size. For facetted interfaces, that typically prevail in vapor and solution growth systems, 
and on which strongly anisotropic growth occurs via atomic layer spreading, the necessary 
conditions for morphological stability are only partly understood. Several workers have proposed an 
isotropic continuum formulation for morphological stability in vapor growth. These models are 
based on the assumption that stability is controlled either by the temperature gradient [25, 26] or by 
capillarity and the dependence of equilibrium concentration on curvature [27]. However, the 
morphologies of vapor-solid [28-32] and solution-solid interfaces [2, 33] are found to be more stable 
than predicted by isotropic models. This suggests that interface kinetics play a key role in 
morphological stability during facetted growth. It has been argued that polycrystalline vapor-solid 
interfaces should exhibit isotropic morphological stability [27]. In general, however, this argument is 
incorrect, for as long as the grain or facet size exceeds the adatom surface diffusion distance of, 
typically, 100 A, anisotropic layer growth will prevail [3, 28]. The stabilizing effect of atomistic 
interfacial kinetics has been incorporated in several anisotropic morphological stability models [34- 
41] Again, probably the most important result obtained from these models is the loss of stability 
beyond* a critical crystal size [34. 41], Such behavior has been observed in experiments [29-33]. Yet 
a complete solution which includes coupling between nutrient transport, the moving boundary and 
anisotropic surface kinetics is still lacking. 

An alternative to continuum descriptions of the growth process is the tracking of individual 
growth units using a Monte Carlo (MC) simulation method [5,42]. MC simulations have been used 
in a number of studies of the dependence of surface roughening and growth rate on supersaturation, 
and temperature [5,43,44]. These studies are primarily focused on microscopic interfacial kinetics, 
and occasionally include surface diffusion [45]. The influence of transport of growth units in the 
bulk nutrient to the interface has received less attention. In the event that interface kinetics govern the 
growth morphology, this simplification is not a severe limitation. However, in reality, nutrient bu 
transport also appears to play a decisive role and often limits the stability of growing faces [29-33, 

46]. 


The recently developed diffusion-limited aggregation (DLA) model of Witten and Sander [47] 
seems particularly appropriate to describe nutrient diffusion in crystal growth. The DLA model 
simulates the growth of an aggregate by considering the random walk of a succession of particles on 


a discretized grid. The DLA model has been adopted by many workers to simulate various types of 
non-equilibrium pattern formation [48, 49]. Typical modifications include the introduction of 
"capillarity" [50, 51], anisotropic sticking probability [52, 53] and noise reduction [54-57]. 
Recently, we have modified the DLA model to include anisotropic interfacial attachment kinetics and 
surface diffusion to study the morphological evolution of crystals growing from a stagnant vapor or 
from a vapor with externally imposed uniform and rotational drift [58-61]. A variety of conditions, 
ranging from kinetics- to transport-controlled growth were examined, and successive transitions 
from compact to open dendritic morphologies were obtained. Diffusion in the nutrient, and increases 
in growth temperature and bulk supersaturation were found to be morphologically destabilizing, 
whereas surface diffusion and anisotropy in attachment kinetics act stabilizing. Specifically, we 
found that a given set of growth parameters defines a critical crystal size beyond which the crystal s 
faces cannot retain their macroscopic flatness, i. e. become morphologically unstable. Surface 
diffusion, on the other hand, continues to play a stabilizing role for micro-facets long after the loss of 
macroscopic stability. In the following section we will summarize our efforts and results. Most 
recently, Saito and Ueta have obtained similar growth shape evolutions from a Monte Carlo model 
with anisotropic attachment and evaporation kinetics features, using a DLA techmque which launches 
many walkers at once [62]. 

MODEL 

Crystal growth consists, in general, in the transport of growth units towards an interface, followed 
by a typically rather complex series of interfacial kinetics processes, before final incorporation into 
the lattice occurs [63]. The rates of the interface kinetics processes of individual growth units will be 
determined by the local configuration of the interface sites that the units happen to visit. Hence, at 
molecular length scales, the crystal symmetry will inevitably be reflected in the anisotropy of the 
attachment kinetics. This underlying symmetry will also be manifested at macroscopic length scales 
unless the atomic roughness of the interface is high enough to allow diffusion in the nutrient to 
control the evolution of the growing shape. Though a complete model of such complex scenarios is 
not practical at this point, we have formulated a Monte Carlo model which retains the essential 
physics of both nutrient transport and interface kinetics, including surface diffusion. 

The gaseous nutrient is assumed to consist of highly dilute growth species A in an inert gas B 
that randomizes the motion of A. The random walk (diffusion) of A in the nutrient is restricted to a 
discrete lattice of spacing b, with equal jump length (representative of a mean free path) , a, in all 
directions, unless a drift is superimposed on the random walk. The grid spacing b represents a 
lattice unit of the crystalline nucleus, which is initially prescribed in the central region of the nutnent, 
and which grows through the addition of random walkers. While we may choose the jump length 
and grid spacing to be the same to minimize computation time, in a real vapor-crystal system a »b. 

To approximate the complex processes following the arrival of a growth unit (random 
walker) at the interface we made the following assumptions. The impingement rate can be obtained, 
based on ideal gas kinetics, from the chemical potential difference between supersaturated bulk vapor 

and average surface site, Aji, in the form [58, 61] 


K + = K eq exp (A(i/kT) , 


( 1 ) 



where K e q is the temperature dependent equilibrium value of K + . Since the evaporation rate is 
sensitive to the local configuration of the site from which a unit is to be dislodged, we cast it, 
following Gilmer and Bennema [45], into the site-dependent form 

K[ = v exp(-Ei/kT) , (2) 

where v is a lattice vibration factor and Ej is the product of the pair interaction (bond) energy of a 
unit with a nearest neighbor and the number of occupied neighbor sites of site i, n]. 

The probability that a growth unit sticks onto a crystal surface is then P] = K + /(K + +Kj). In 
order to relate Keq and v we assume local equilibrium and thus equate impingement and evaporation 
rates at equilibrium and obtain [58] 

p _ 

1+7 p (n °' n i ) t (3) 

where no, the number of nearest neighbors in a kink site, is 2 for a square lattice and 3 for a 
triangular lattice, and v = exp (A|i/kT) and [3 = exp (-<))/kT). 

If the growth unit does not stick onto the site i, it will jump to one of its unoccupied neighbor 
sites on the surface or in the nutrient. The jump probability from site i to a neighboring unoccupied 
sitej is assumed to be [58, 61] 



(4) 


where c' is the number of unoccupied nearest neighbor sites of site i. Clearly, a larger nj results in a 
higher probability that a molecule will jump to site j on the interface. 

The combination of equations 3 and 4 contains the essential physics of the surface kinetics 
for a crystal growing from a vapor. It should be realized that the crystal anisotropy has already been 
taken into account implicitly in equations 3 and 4, since both nj and nj are anisotropic factors. 


The strength ci surface diffusion (described by equation 4) is characterized by the diffusion 
time or diffusion length on the surface. An exact calculation^ surface diffusion length is difficult. In 

our model we estimate the average surface diffusion length A, s according to 


X s ~ b , (5) 

where n t is the average (residence) "time" (actually the number of MC steps) for a random walker to 
diffuse on the interface. 




Figure 1. Flowchart of simulation steps. 


SIMULATION PROCEDURE 

A summary of the steps involved in the MC simulation is presented in the flowchart of figure 
1. A growth unit is released from a circular source that surrounds the prescribed crystal nucleus. 
(The choice of source diameter is discussed in [58, 61]). The random walk is continued until the 
walker encounters the interface. The number of neighbors of the impingement site is then 
determined and the sticking probability Pi calculated according to equation 3. A random number R, 
0 <R < 1, is generated. If R < Pi, a registration to site i is made, and then that walker is retired. 
When an interfacial site has been registered N times it is considered to be occupied. A larger number 
of registrations N has the effect of reducing the noise inherent in a simulation of this kind [54-58]. If 









R > Pi, the walker will leave site i without registration and jump to one of its unoccupied 
neighboring sites. The jump probability which determines the new location on the interface or in the 
nutrient is calculated according to equation 4 and compared to a random number R’. As the walker 
jumps to a new unoccupied site i\ the above procedure is repeated. For typical values of the 
interaction parameters the probability that the particle returns to the nutrient is much less than that for 
continued residence on the surface. Successive walkers are released from the outer source and they 
proceed to undergo random motion until they either escape across the source boundary or make a 
registration at an interfacial site. The simulation is continued until the crystal reaches a specified size 
M. Note that, in contrast to real systems, once a site is considered occupied, its particle cannot leave 
the site, i.e. evaporation is only permitted in transit. 

The simulation was implemented on a CRAY X-MP/24. Typical CPU times used for each 
simulation were between 2,000-7,000 sec. To show the sequence of growth patterns, we divided 
the total growth particles into consecutive groups and outlined the morphology resulting after 

completion of each group. 

results and discussion 

Growth morphology and nutrient co ncentration 
distribution at mean free p ath of one lattice_unil 

First we present selected results obtained with a = b. In this case the simulation was terminated 
when the growth pattern had gained 3xl0 3 particles, which were consecutively divided into tlvee 
groups of 1,000 each. Figure 2 shows the effect of growth temperature and/or bond strength on 



Figure 2. Effect of bond strength and temperature on growth patterns of a crystal (triangular lattice) 
at a fixed normalized supersaturation (Ap/kT = 0.69). (a) <(>/kT = 3.91; (b) <j)/kT = 2.30, (c) <(>/kT - 
0.69; (d) <t>/kT = 0.36. 



growth morphology at a fixed normalized supersaturation (Ap/kT = 0.69 or y = 2.0). At low 
temperatures the growing "crystals" retain a compact, facetted form up to their terminal size (figure 
2a) despite the appreciable supersaturation. With increasing temperature or decreasing <{>/kT, the 
even smaller growth forms, i.e. after 1000-2000 growth units, show "microscopic" depressions (of 
a few lattice constants) in the center of the facets. Later growth under these or even lower <}>/kT 
conditions, leads to the evolution of "macroscopic" depressions (figure 2b). On further increase of 
the temperature, this transition occurs at even smaller sizes and dendritic growth evolves (figure 2c) 
which subsequently exhibits extensive side-branching (figure 2d). 

Figure 2 provides considerable insight into the morphologically stabilizing effect of 
anisotropic surface kinetics. One sees that the formation of depressions always sets in at facet (face) 
centers. This, as has been shown in detail in [60], is associated with a local decrease in both the 
nutrient concentration and its gradient. However, at the comers, which are more readily reached by 
building blocks, both the concentration and its gradient remain high. Furthermore, as can be seen 
from the sequence of contours of the lower temperature parts of figure 2, the formation of a 
microscopic depression is associated with the exposure of a higher step density in that region. 
Hence, both local sticking probability (equation 3) and jumping probability from neighboring sites 
(equation 4) increase. Consequently, not only is the leaner supply of growth units in the central 
region better utilized, but the surface diffusion flux is also increased. Both aid the crystal in retainmg 
its (macroscopically) facetted morphology. This stabilizing effect reaches a limit as, with increasing 
depression, the growth step density becomes sufficient to accommodate the attachment of most 
growth units that arrive at the depleted face centers. Further increase in the crystal size is then 
accompanied by the loss of morphological stability. 

This interplay between destabilizing bulk diffusion and stabilizing interface kinetics has been 
invoked before in Chernov's theory of growth-shape preservation for facetted crystals [34] and its 
adaption by Kuroda et al. [35]. 

The above simulation also shows, for the first time, the destabilizing effect of an increase in 
the growth temperature, that can be deduced from the transition from facetted to non-facetted shapes 
in figure 2. This transition is loosely speaking analogous to thermal roughening; for a detailed 
discussion see [64-66]. 

We have also investigated the effect of supersaturation on growth morphology [58], which 
enters in our model through the dependence of the sticking probability Pi (equation 3) on the 
impingement rate K + (equation 1). We found that an increase in supersaturation also decreases the 
critical crystal size beyond which the macroscopic face stability is lost. This is analogous to kinetic 
roughening observed in earlier MC modelling [67-69] and in experimental studies [69-73]. The 
dependence of the crystal morphology on normalized supersaturation and bond strength/temperature 
is tightly coupled. All our results for a triangular lattice are summarized in figure 3 [58]. Depending 
on the value of <|>/kT and Ap/kT, the growing crystal can acquire one of the following morphologies: 
compact facetted (region I), compact branched (region II), and open branched and sidebranched 
(region III). In general, the supersaturation above which the compact facetted form is no longer 
stable decreases with decreasing bond strength or increasing temperature. 




Figure 3. Crystal-growth morphology as a collective effect of supersaturation and bond strength (or 
temperature) for triangular lattice at mean free path of one lattice unit. In region I crystals acquire 
compact-facetted forms; in region II compact branched with six-fold symmetry, in region T I 
dendritic with multiple sidebranches. Symbols: Monte Carlo results judged as boundary cases . 
Solid line: Temkin's boundary between facetted and nonfacetted growth for a simple cubic lattice 
[74]. For further details see [58]. 

Changes in growth morphology with in c rease in mean free path 


In this section we increase the mean free path and present growth patterns arising on a square 
lattice after addition of four groups of 2,500 particles each. Figure 4 shows the effect on 
morphology of bond strength or temperature at a constant normalized supersaturation, and a mean 
free path of a = 10b. As compared to results obtained on a square lattice for a=b [60] the critical size 
of a facet has increased about tenfold. The stabilizing effect of an increase in the mean free path is 
further illustrated in figure 5, which is for the same fixed values of growth temperature and 
supersaturation as figure 4d. One sees that the growth morphology becomes increasingly compact 
and facetted as a is increased in steps from 5 to 100b. Such morphological changes have been 
observed experimentally [29-33]: growing facetted crystals became centrally depressed as the mean 
free path was decreased by increasing the partial pressure of the inert gas. We have found that the 
critical crystal size scales linearly with the mean free path in the parameter range considered [60]. 
Extrapolation of the critical size to more realistic conditions [60] gave reasonable agreement with 
critical sizes found experimentally [31]. 

Critical inspection of figures 4 and 5 shows that even under conditions in which macroscopic 
facet stability is lost, the contours possess microscopic facets, i.e. flat regions of a much shorter 
length scale which, as shown in [60, 61], corresponds to the average surface diffusion length. 








Figure 4. Effect of bond strength and temperature on growth patterns at a large mean free path (a = 
10 lattice units) at a fixed super saturation (A|i/kT = 0.69). (a) <J>/kT = 4.60; (b) <J)/kT = 3.91; (c) 
4>/kT = 2.30; (d) <J>/kT = 0.69. Contours corresponds to addition of 2,500 particles each. 


Changes in crystal morphology in t he presence of a drift 

With an externally imposed drift, the individual jumps of a random walker will be biased, see 
[59, 60]. Figure 6 shows the effect of a uniform drift (i.e. with an additional unidirectional 
displacement of d = 2b at each random jump) on growth morphology for the same conditions as 
figure 4. One sees that the drift causes asymmetric growth. Sites facing the drift acquire building 
blocks more readily than those on the lee side. Although the drift affects all cases, the response of 
the growth morphology increases with decreasing bond strength or increasing temperature (figure 
6a-d). In addition to modifying the growth rate, the drift is also seen to stabilize the front face in 
comparison to the corresponding "no-drift cases of figure 4. 

In addition to unidirectional drifts, we have also investigated the consequences of rotational 
drifts [60]. In particular, rotational drifts with periodic reversal of direction were found to stabilize 
facetted morphologies. This is in agreement with observed growth morphologies of crystals subject 
to frequent rotation reversals [60]. 





Figure 5. Effect of mean free path on 
growth patterns at A(i/kT = 0.69 and 
<j)/kT = 0.69. (a) a = 5b; (b) a = 25b; (c) 
a = 50b; (d) a= 100b. Contours 

correspond to addition of 2,500 particles 
each. 


(C) 


(d) 



Figure 6. Effect of drift with d = 2b on 
growth morphology. Growth 
conditions otherwise the same as in 
figure 4, i.e., A|i/kT = 0.69. (a) 0/kT 
= 4.60; (b) <(>/kT = 3.91; (c) <J)/kT = 
2.30; (d) <>/kT = 0.69. The arrow 
indicates the direction of the drift. 


(a) 
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